I am trying to make consensus gene sequences for multiple gene queries from a sample of forward and reverse paired end genomic reads. As of now I am just testing on 5 - 7 reference gene sequences, but plan to use it for a MUCH larger set and with multiple samples.
For each gene query, the forward reads are blasted with "blastall" against the refence gene sequence (gene specific single fasta record) and the reads with highly significant alignments are written to a file (geneA.FORWARD.fastq) and their mate pairs from the reverse reads are written to geneA.REVERSE.fastq. The reverse reads are also blasted and significant hits are written to geneA.FORWARD.fastq with the mate pairs from forward reads written to geneA.REVERSE.fastq. Now one file has all the hits from f/r reads and the other has all the mate pairs.
I then use bwa (0.7), samtools (1.1; tried 1.2 as well), and bcftools (1.2; tried 1.1 as well) to index the reference file (contains fasta records for all reference genes used with each record being identical to the one used for the blast step), create an alignment with the reference and reads, and finally a consensus sequence.
The commands used are:
On some of the genes it works just fine and I get my consensus sequence for that gene. However, on others, I would get the error "Segmentation fault" on the very last step (bcftools consensus).
After backtracking for the issue I changed the samtools mpileup step to
Here I noticed that for the genes that worked, the stdout was the vcf header followed by data, but for the genes that didn't it was just the vcf header with no data below. So it seems the issue comes before or during this step from what I was originally doing. Another head scratcher is that for the genes that didn't work, if I take out the mate pairs for the alignment step and do
or
it works just fine.
Any ideas on what I can look for with the reads of the genes that didn't work? All the queries seem to have good, significant blast hits being used and everything is being done the same way for working queries and non working queries. All the gene sequences are very similar and the files with the reads with significant hits all seem fine. After the bwa mem step, the sam files created seem fine, but I am not really sure what I can look for that would show what is different from the ones that did work and the ones that didn't.
For each gene query, the forward reads are blasted with "blastall" against the refence gene sequence (gene specific single fasta record) and the reads with highly significant alignments are written to a file (geneA.FORWARD.fastq) and their mate pairs from the reverse reads are written to geneA.REVERSE.fastq. The reverse reads are also blasted and significant hits are written to geneA.FORWARD.fastq with the mate pairs from forward reads written to geneA.REVERSE.fastq. Now one file has all the hits from f/r reads and the other has all the mate pairs.
I then use bwa (0.7), samtools (1.1; tried 1.2 as well), and bcftools (1.2; tried 1.1 as well) to index the reference file (contains fasta records for all reference genes used with each record being identical to the one used for the blast step), create an alignment with the reference and reads, and finally a consensus sequence.
The commands used are:
Code:
bwa index ref.fa bwa mem -t 4 ref.fa geneA.FORWARD.fastq geneA.REVERSE.fastq > geneA.sam samtools view -S -b -F4 geneA.sam > geneA.bam samtools sort geneA.bam Sorted.geneA samtools index Sorted.geneA.bam samtools mpileup -vf ref.fa Sorted.geneA.bam > Sorted.geneA.vcf.gz bcftools index Sorted.geneA.vcf.gz bcftools view -r geneA:start-end -o Extracted.geneA.vcf.gz -O z Sorted.geneA.vcf.gz bcftools index Extracted.geneA.vcf.gz bcftools call -c -r geneA:start-end -A -O z -o Call.geneA.vcf.gz Extracted.geneA.vcf.gz bcftools index Call.geneA.vcf.gz samtools faidx ref.fa geneA:start-end | bcftools consensus Call.geneA.vcf.gz > cns.geneA.fa
After backtracking for the issue I changed the samtools mpileup step to
Code:
samtools mpileup -uf ref.fa Sorted.geneA.bam | bcftools call -c
Code:
bwa mem -t 4 ref.fa geneA.FORWARD.fastq > geneA.sam
Code:
bwa mem -t 4 ref.fa geneA.REVERSE.fastq > geneA.sam
Any ideas on what I can look for with the reads of the genes that didn't work? All the queries seem to have good, significant blast hits being used and everything is being done the same way for working queries and non working queries. All the gene sequences are very similar and the files with the reads with significant hits all seem fine. After the bwa mem step, the sam files created seem fine, but I am not really sure what I can look for that would show what is different from the ones that did work and the ones that didn't.