#Amomum


#Adaptor and quality trimming using Trimmomatic
java -jar trimmomatic-0.38.jar PE -phred33 Amomum-subulatum_S51_L001_R1_001_22000.fastq.gz Amomum-subulatum_S51_L001_R2_001_22000.fastq.gz Amomum-1P.fastq.gz Amomum-1U.fastq.gz Amomum-2P.fastq.gz Amomum-2U.fastq.gz ILLUMINACLIP:adapters\TruSeq3-PE.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:5:20 MINLEN:36


#Duplicate removal using fastuniq
#unzip
gzip -d Amomum-1P.fastq.gz
gzip -d Amomum-2P.fastq.gz
#prepare list of sequences
ls Amomum-?P.fastq > fastqlist.txt
#run duplicate removal
fastuniq -i fastqlist.txt -t q -o Amomum-1P_nodup.fastq -p Amomum-2P_nodup.fastq

#count number of reads in FASTQ file (divide the resulting numbers by 4) (wc -l = counts number of lines)
#original file must be unzipped first
gzip –d Amomum-subulatum_S51_L001_R2_001_22000.fastq.gz
#before quality trimming
cat Amomum-subulatum_S51_L001_R2_001_22000.fastq | wc -l
#after quality trimming
cat Amomum-1P.fastq | wc -l
#after duplicate removal
cat Amomum-2P_nodup.fastq | wc -l

#gzip the deduplicated samples (again, '?' allows to run the command just once)
gzip Amomum-?P_nodup.fastq


#Read mapping with BWA
#index the reference
bwa index curcuma_HybSeqProbes_test_with400Ns_beginend.fas
#map reads to the reference
bwa mem curcuma_HybSeqProbes_test_with400Ns_beginend.fas Amomum-1P_nodup.fastq.gz Amomum-2P_nodup.fastq.gz > Amomum.sam
#convert to BAM
samtools view –bS Amomum.sam > Amomum.bam
#count number of all reads in a BAM file
#cut - takes 1st column/field, i.e., read name
#sort - alphabetically sort lines (necessary to do before applying the next command)
#uniq - report unique names only, i.e. remove duplicated names (reads can mapped to multiple locations)
#wc - 'word count', counts number of lines with the parameter '-l'
samtools view Amomum.bam | cut -f1 | sort | uniq | wc -l
#count number of alignments in BAM file (reads mapped to multiple locations count multiple times)
samtools view -c Amomum.bam
#count number of mapped alignments (flag 'unmapped' absent, i.e. the 3rd bit - binary = 4, hexadecimally = 0x40 - is not set - '-F')
samtools view -F 0x04 -c Amomum.bam
#count number of unmapped alignments (flag 'unmapped' present - '-f')
samtools view -f 0x04 -c Amomum.bam
#count number of supplementary alignments (flag 'supplementary alignment' present, i.e. the 12th bit - binary = 2048, hexadecimally = 0x800 - is set)
samtools view -f 0x800 -c Amomum.bam
#count number of mapped reads (alignments not flagged as unmapped (0x4) or supplementary (0x800))
samtools view -F 0x804 -c Amomum.bam
#sort the BAM file
samtools sort Amomum.bam Amomum_sorted
#index the sorted BAM file
samtools index Amomum_sorted.bam
#extract only mapped reads (reads must be flagged as mapped in proper pair (0x2))
samtools view -u -f 0x2 Amomum_sorted.bam > Amomum_mapped.bam
#sort BAM file with extracted reads by the read name
samtools sort -n Amomum_mapped.bam Amomum_mapped.sort
#extract mapped reads as pair-end FASTQ using bedtools
bedtools bamtofastq -i Amomum_mapped.sort.bam -fq Amomum_read1_mapped.fastq -fq2 Amomum_read2_mapped.fastq
#gzip the resulting FASTQ files
gzip Amomum_read?_mapped.fastq


#Variant calling using SAMtools/BCFtools
#create 'pileup' for the sorted BAM file (u = uncompressed output)
samtools mpileup -uf curcuma_HybSeqProbes_test_with400Ns_beginend.fas Amomum_sorted.bam > Amomum.pileup
#call variants/genotypes using bcftools (v = variant sites only, c = SNP calling, g = call genotypes) [this is the command for old version, the new version uses 'bcftools view' and difeerent parameters]
bcftools view -vcg Amomum.pileup > Amomum.vcf
#how many variants are in VCF file (grep -v = prints all lines not containing a string, "^#" - line starting with #, wc -l = count number of lines)
grep -v "^#" Amomum.vcf | wc -l
#how many homozygotes (1/1) and heterozygote (0/1) are present?
grep "1/1" Amomum.vcf | wc -l
grep "0/1" Amomum.vcf | wc -l
#filter only variants with quality (QUAL) of at least 36 (awk deals with columns, i.e. if the value in 6th column ($6) is >= 36 it prints the whole line - $0)
awk "{if ($6 >= 36) print $0}" Amomum.vcf | grep "0/1" | wc -l
#filter only variants covered by at least 15 reads - DP in the INFO column is read depth (awk -F"=|;" - sets the column separator to both '=' and ';' - this makes the DP value 2nd column
awk -F"=|;" "{if ($2 >= 15) print $0}" Amomum.vcf | grep "0/1" | wc -l
