---
title: "transfer"
format: html
editor: visual
---

# Code for running subset of files on Raven
Preliminary work on for running full alignment on 1 _Prochlorococcus_ genome with 1 G3 diel sample

## Transfer files from armbrust lab server
Previously, I set up an SSH key pair in order to log in to my server from this server. These files have already been trimmed and Qc'd.

```{bash}
ls
```


```{bash}
## grab a subset of files from grazer (1 for now)
# only want diel experiments and paired data

# grab only a subset of files and save in txt file
ssh klqi@grazer "ls /mnt/nfs/projects/armbrust-metat/gradients3/g3_station_ns_metat/qc_data/G3.DIEL*.fastq.gz | grep -v '.unpaired.fastq.gz'" | head -n 6 > file_list.txt

# recursively grab files from grazer
while read file; do
    rsync -av klqi@grazer:"$file" /home/shared/16TB_HDD_01/fish546/kathy-pico/data/trimmed/
done < file_list.txt
```

# Align fastq files to reference genomes

https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000015665.1/ Try one first, then use all \~600 from stephen's paper https://www.science.org/doi/10.1126/sciadv.adr4310, https://www.nature.com/articles/sdata2018154#Sec2

```{bash}
# grab reference genome from NCBI
cd ~/github/project/alignment
wget -L https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/015/665/GCF_000015665.1_ASM1566v1/GCF_000015665.1_ASM1566v1_genomic.fna.gz
# change the name
mv GCF_000015665.1_ASM1566v1_genomic.fna.gz mit_9515_reference.fna.gz
# unzip
gunzip mit_9515_reference.fna.gz
```

```{bash}
## try one first (bowtie)
# build bowtie index
cd ../alignment
# mkdir output
cd output
bowtie-build mit_9515_reference.fna mit_9515

# grab files from data folder
datadir="/home/shared/16TB_HDD_01/fish546/kathy-pico/data/trimmed/"
# create paths from first 2 lines of file
read1=$datadir"G3.DIEL.NS.S4C11.15m.0_2um.A.fw.fastq.gz"
read2=$datadir"G3.DIEL.NS.S4C11.15m.0_2um.A.rv.fastq.gz"
# align with bowtie
bowtie -q mit_9515 -1 "$read1" -2 "$read2" -S test_alignment.sam
```

Output with bowtie and MIT 9515:

```         
# reads processed: 67980220
# reads with at least one reported alignment: 9122 (0.01%)
# reads that failed to align: 67971098 (99.99%)
Reported 9122 paired-end alignments
```

Maybe that wasn't the correct reference genome, I should probably try a different ecotype - Actually, Zinka said this was not bad since I just used 1 genome and got \>9000 reads that aligned

# process .sam file

```{bash}
# visualize the aligments

```

Maybe it would be better to include the annotations with bowtie

# Now try with bowtie2

```{bash}
# export path to bowtie2 to run scripts
bowtie2dir="/home/shared/bowtie2-2.4.4-linux-x86_64/"
cd $bowtie2dir
# build index with bowtie2
mkdir ~/github/project/alignment/output2
# store files in output2 folder
./bowtie2-build ~/github/project/alignment/mit_9515_reference.fna \
~/github/project/alignment/output2/mit_9515_reference2

# grab files from data folder
datadir="/home/shared/16TB_HDD_01/fish546/kathy-pico/data/trimmed/"
# create paths from first 2 lines of file
read1=$datadir"G3.DIEL.NS.S4C11.15m.0_2um.A.fw.fastq.gz"
read2=$datadir"G3.DIEL.NS.S4C11.15m.0_2um.A.rv.fastq.gz"
# align with bowtie2
./bowtie2 -x ~/github/project/alignment/output2/mit_9515_reference2 \
-1 $read1 \
-2 $read2 \
-S ~/github/project/alignment/output2/test_alignment2.sam
```

Output with Bowtie2 and MIT9515:

```         
67980220 reads; of these:
  67980220 (100.00%) were paired; of these:
    66197024 (97.38%) aligned concordantly 0 times
    1783114 (2.62%) aligned concordantly exactly 1 time
    82 (0.00%) aligned concordantly >1 times
    ----
    66197024 pairs aligned concordantly 0 times; of these:
      22993 (0.03%) aligned discordantly 1 time
    ----
    66174031 pairs aligned 0 times concordantly or discordantly; of these:
      132348062 mates make up the pairs; of these:
        130429252 (98.55%) aligned 0 times
        1917243 (1.45%) aligned exactly 1 time
        1567 (0.00%) aligned >1 times
4.07% overall alignment rate
```

# look at .sam file

```{bash}
cd ~/github/project/alignment/output2
# get some statistics
samtools flagstat test_alignment2.sam

# convert to BAM file- i think this is better but i don't know why
samtools view -bh test_alignment2.sam > test_alignment2.bam
samtools sort test_alignment2.bam -o test_alignment2_sorted.bam

# grab annotation file from ncbi
wget -L https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/015/665/GCF_000015665.1_ASM1566v1/GCF_000015665.1_ASM1566v1_genomic.gff.gz -O mit9515_annotation.gff.gz
# unzip
gunzip mit9515_annotation.gff.gz

# get counts (featureCounts?)
```

Statistics Output:

```         
135960440 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
5531188 + 0 mapped (4.07% : N/A)
135960440 + 0 paired in sequencing
67980220 + 0 read1
67980220 + 0 read2
3566392 + 0 properly paired (2.62% : N/A)
3612380 + 0 with itself and mate mapped
1918808 + 0 singletons (1.41% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
```


## Extract only predicted gene counts for all pro genomes, instead of full genome