Variant Calling for Cafeteria roenbergensis strain RCC970

Download data

Download the reference genome - Cafeteria roenbergensis strain BVI as archive containing the genome in FASTA format from NCBI https://www.ncbi.nlm.nih.gov/assembly/GCA_008330645.1/

wget -O BVI.tar https://www.ncbi.nlm.nih.gov/projects/r_gencoll/ftp_service/nph-gc-ftp-service.cgi/?HistoryId=MCID_60dfdfa5ccd4273699521e5b&QueryKey=2&ReleaseType=GenBank&FileType=GENOME_FASTA&Flat=true
# uncompress file
tar -xf BVI.tar
# place uncompressed FASTA file at top directory for easy access
gzip -cd *ncbi-genomes-2021-07-03/GCA_008330645.1_CrBVI_genomic.fna.gz > BVI.fna
# remove temporary files
rm -r ncbi-genomes* 

Download paired-end sequencing data - Cafeteria roenbergensis strain RCC970 - in FASTQ format from NCBI - https://www.ncbi.nlm.nih.gov/sra/SRX6482057[accn].

fastq-dump --split-e SRR9724623
# Read 6619536 spots for SRR9724623
# Written 6619536 spots for SRR9724623

Check read quality

Check read quality, trim & filter with `fastp`. Takes about 1-2 minutes. The results are summarized in fastp.html.

fastp -w 60 --cut_right -i SRR9724623_1.fastq -I SRR9724623_2.fastq -o RCC_1.fq -O RCC_2.fq

Map reads

Map each read pair to the optimal location on the reference assembly with `bwa mem` and save in `.bam` format. Takes about 2-10 minutes depending on machine.

bwa index BVI.fna
bwa mem -t 60 BVI.fna RCC_1.fq RCC_2.fq > BVI-RCC.sam
samtools view -ub BVI-RCC.sam | samtools sort -@ 60 -m 15G -o BVI-RCC.bam /dev/fd/0

We are foregoing PCR-duplicates marking and removal because a) we used a protocol with highly reduced PCR-cycles and b) Ebbert et al. (2016) BMC Bioinformatics showing that PCR-dedup has no effect of variant calls.

Call & filter variants

freebayes -f BVI.fna BVI-RCC.bam > BVI-RCC.vcf
bcftools filter --include 'QUAL>20 & INFO/DP>30 & INFO/DP < 150' BVI-RCC.vcf >BVI-RCC-filtered.vcf