#Ensete

#Download FASTQ files from SRA
fastq-dump --split-3 -gzip SRR7058210


#Adaptor and quality trimming using Trimmomatic
java -jar trimmomatic-0.38.jar PE -phred33 Ensete-superbum_S375_L001_R1_001.fastq.gz Ensete-superbum_S375_L001_R2_001.fastq.gz Ensete-1P.fastq.gz Ensete-1U.fastq.gz Ensete-2P.fastq.gz Ensete-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 Ensete-1P.fastq.gz
gzip -d Ensete-2P.fastq.gz
#prepare list of sequences
ls Ensete-?P.fastq > fastqlist.txt
#run duplicate removal
fastuniq -i fastqlist.txt -t q -o Ensete-1P_nodup.fastq -p Ensete-2P_nodup.fastq

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

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


#Read mapping with BWA
#index the reference
bwa index HF677508_Musa_acuminata_cpDNA.fas
#map reads to the reference
bwa mem HF677508_Musa_acuminata_cpDNA.fas Ensete-1P_nodup.fastq.gz Ensete-2P_nodup.fastq.gz > Ensete.sam
#convert to BAM
samtools view -bS Ensete.sam
#count number of all reads in a BAM file 
samtools view Ensete.bam | cut -f1 | sort | uniq | wc -l #not working on Windows
#samtools view Ensete.bam | cut -f1 > enstemp #only for Windows
#sort < enstemp | uniq | wc -l #only for Windows
#count number of alignments in BAM file (reads mapped to multiple locations count multiple times)
samtools view -c Ensete.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 Ensete.bam
#count number of unmapped alignments (flag 'unmapped' present - '-f')
samtools view -f 0x04 -c Ensete.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 Ensete.bam
#count number of mapped reads (alignments not flagged as unmapped (0x4) or supplementary (0x800))
samtools view -F 0x804 -c Ensete.bam
#sort the BAM file
samtools sort Ensete.bam -o Ensete_sorted.bam
#index the sorted BAM file
samtools index Ensete_sorted.bam
#extract only mapped reads (reads must be flagged as mapped in proper pair (0x2))
samtools view -u -f 0x2 Ensete_sorted.bam > Ensete_mapped.bam
#sort BAM file with extracted reads by the read name
samtools sort -n Ensete_mapped.bam -o Ensete_mapped.sort.bam
#extract mapped reads as pair-end FASTQ using bedtools
bedtools bamtofastq -i Ensete_mapped.sort.bam -fq Ensete_read1_mapped.fastq -fq2 Ensete_read2_mapped.fastq
#gzip the resulting FASTQ files
gzip Ensete_read?_mapped.fastq


#De-novo assembly using Velvet
#run velveth with k-mer 101 (k101 - name of the output directory, 101 - k-mer length). K-mer has to be an odd number
velveth k101 101 -fastq -shortPaired -separate Ensete_read1_mapped.fastq.gz Ensete_read2_mapped.fastq.gz
#run velvetg on the previous results
velvetg k101 -cov_cutoff auto -min_contig_lgth 1000 -exp_cov auto
