Scripts used for assembly of RNASE1 genes.
Janiak MC, Burrell AS, Orkin JD & Disotell TR. (2019). Duplication and convergent evolution of the pancreatic ribonuclease gene (RNASE1) in folivorous non-colobine primates, the howler monkeys (Alouatta spp.). Scientific Reports 9, 20366. https://doi.org/10.1038/s41598-019-56941-7
We assembled the RNASE1 gene region for four howler monkey species, from whole genome short reads generated with the 10X Genomics system and from short reads generated by sequencing amplicons on a MiSeq.
To do a targeted re-assembly of a region of interest with convential tools, it is necessary to remove the 10X linked barcodes from the reads, which we did with the script process_10xReads.py (from https://github.com/ucdavis-bioinformatics/proc10xG).
python ~/scripts/process_10xReads.py -a -o ./trimmed_fastq/ -1 ./fastq/HowlerMonkey_S56_L007_R1_001.fastq.gz -2 ./fastq/HowlerMonkey_S56_L007_R2_001.fastq.gz
Below is the code for running each step of the assembly process for one of our samples. In reality, these were sent out as jobs and run in parallel on the Amarel cluster at Rutgers University and on the Helix cluster at the University of Calgary. Job scripts are in the scripts directory.
10X short reads with barcodes removed and the demultiplexed MiSeq reads were mapped to the human RNASE1 region with bwa mem and sorted and indexed with samtools (https://github.com/lh3/bwa)
bwa mem -M ~/Alouatta_mapping/human_RNASE1_fit.fa ~/Alouatta_mapping/MiSeq_7March19/Apal_S1_L001_R1_001.fastq.gz ~/Alouatta_mapping/MiSeq_7March19/Apal_S1_L001_R2_001.fastq.gz -t4 | samtools view -@4 -Sbh - | samtools sort -@4 -T ~/Alouatta_mapping/mapping/bwa_out/Apal_S1.sort.tmp -O bam -o - > ~/Alouatta_mapping/mapping/bwa_out/Apal_S1.sort.bam; samtools index ~/Alouatta_mapping/mapping/bwa_out/Apal_S1.sort.bam
After this initial mapping to the human RNASE reference, we repeated the mapping step using the consensus sequence from the bam files output from the previous mapping step for each species. Consensus sequences were exported from IGV.
bwa mem -M ~/Alouatta_mapping/mapping/consensus_sequences/Apal_S1.consensus.fa ~/Alouatta_mapping/MiSeq_7March19/Apal_S1_L001_R1_001.fastq.gz ~/Alouatta_mapping/MiSeq_7March19/Apal_S1_L001_R2_001.fastq.gz -t4 | samtools view -@4 -Sbh - | samtools sort -@4 -T ~/Alouatta_mapping/mapping/bwa_out/Apal_S1.sort.tmp -O bam -o - > ~/Alouatta_mapping/mapping/bwa_out/Apal_S1.consensus.sort.bam; samtools index ~/Alouatta_mapping/mapping/bwa_out/Apal_S1.consensus.sort.bam
To prepare the bam files for variant calling with freebayes, we did the following:
#Add read group
~/bamtools/bamaddrg/bamaddrg -b Apal_S1.consensus.sort.bam -s Apal_S1 > Apal_S1.consensus.RG.bam
#Sort
samtools sort -n -o Apal_S1.consensus.RG.sort.bam Apal_S1.consensus.RG.bam
#Fix paired end reads
samtools fixmate -m Apal_S1.consensus.RG.sort.bam Apal_S1.consensus.RG.fixmate.bam
#Sort again
samtools sort -o Apal_S1.consensus.RG.fixmate.sort.bam Apal_S1.consensus.RG.fixmate.bam
#Mark duplicates
samtools markdup Apal_S1.consensus.RG.fixmate.sort.bam Apal_S1.consensus.RG.markdup.bam
To call variants with freebayes, we first called SNPs and filtered the output by quality score:
~/bin/freebayes -f ~/Alouatta_mapping/variant_calling/freebayes/human_RNASE1_fit.fa ~/Alouatta_mapping/mapping/bwa_out/AddedRG/Apal_S1.RG.markdup.bam --report-all-haplotype-alleles --haplotype-length 0 > ~/Alouatta_mapping/variant_calling/freebayes/Apal_S1.var.vcf
~/freebayes/vcflib/bin/vcffilter -f "QUAL > 20" Apal_S1.consensus.var.vcf > Apal_S1.consensus.var.filtered.vcf
Filtered vcf files were prepped to be used as input for haplotype calling by bgzipping and indexing:
bgzip Apal_S1.consensus.var.filtered.vcf > Apal_S1.consensus.var.filtered.vcf.gz
tabix -p vcf Apal_S1.consensus.var.filtered.vcf.gz
The reported SNPs were used as input to call longer variants for the coding region:
~/bin/freebayes -f ~/Alouatta_mapping/mapping/consensus_sequences/Apal_S1.consensus.fa --region Apal_S1_RNASE1_consensus:690-1155 --haplotype-length 50 --report-all-haplotype-alleles --haplotype-basis-alleles ~/Alouatta_mapping/variant_calling/freebayes/Apal_S1.consensus.var.filtered.vcf.gz ~/Alouatta_mapping/mapping/bwa_out/AddedRG/Apal_S1.consensus.RG.markdup.bam > ~/Alouatta_mapping/variant_calling/freebayes/haplotype_calling/Apal_S1.hap.vcf
Average coverage was calculated with:
samtools depth -a Apal_S1.consensus.RG.markdup.bam | awk '{sum+=$3} END { print "Average = ",sum/NR}'