mapping on to custom reference sequence

[PratikC]

Hi,
I have sequence read data of 3 particular genes of Human.
I know their location on reference genome and so I retrieved their sequences and made another fasta file with regular header (>seq_name) and following concatenated three sequence without anything else.

My raw read file is of 18GB (paired end,Illumina). When I map them to custom reference (using BWA-SAMtools) I get Bam file of 8GB (Forward-Reverse merged). These all steps takes few minutes to few hours on my workstation. My BAM file index is only 36KB. When I input this BAM file in 'samtools pileup -vcf', it keeps running more than 3 days and the output file size is around 1.1GB. Its still running and I dont understand why it takes so much time? Another matter is when I checked output file after first day, it was 900MB, then on second day it was 1 GB and on third day its 1.1GB. I already killed this process before 7-8 days (because of same behavior) and restarted it. Please help me to understand why it is happening so? (there is no any error message on any of the steps)

Thanking you in advance...

[aaronh]

First off, you should most likely be mapping to the whole genome. Even if your data is from something like exon capture, you need to know how likely the mapping is to each position in the whole genome. Second, don't use pileup as this is unsupported now. Use mpileup.

[colindaven]

You have that much sequence of only a few genes ? That seems like a bizarre form of overkill!

I would also align to the whole genome and use mpileup as the previous poster said - there are good instructions on the samtools sourceforge website.

[PratikC]

Sorry friends for late reply. Thnx aaronh and colindaven. Yes my seq read is of 18GB and only 3 genes are there, but keep in mind its pooled sequencing so we got around 25 samples. Now coming to the main question, what I doubt is can pileup -vdf take so much time when rest of the steps are getting done withing few minutes to few hours. Does pileup take most of the ngs data processing time?

Now as you suggested, I got newer samtool and set up in my workstation. Will keep you posted when done with that one. Maybe mpileup will be more efficient to do the thing...

[Heisman]

Pileup should not take that much time. Since you are going to try mpileup, as it is running you can check to see how far it has progressed. Assuming you sort your bam file prior to running it you'll be able to guess how long it will take.

If it's still acting weirdly, write back. I've set up custom reference sequences for this purpose (although I use Novoalign, not BWA, to align) and made it work fine.

[PratikC]

Pileup should not take that much time. Since you are going to try mpileup, as it is running you can check to see how far it has progressed. Assuming you sort your bam file prior to running it you'll be able to guess how long it will take.

If it's still acting weirdly, write back. I've set up custom reference sequences for this purpose (although I use Novoalign, not BWA, to align) and made it work fine.

Hi Heisman, I already started mpileup and hopes it will be completed soon. Meanwhile, if possible, can you plz give me some more details of how you set up custom reference sequences and then how you did SNP annotation? You can also mail me if you want to.
Thanks.

[Heisman]

I actually haven't done this for awhile so hopefully I'm not forgetting anything.

1. Take all of your targeted regions, extend each region by X base pairs, where X is determined on your method of generating data. Put all of this in fasta format.

2. When you have your fasta file, you will want to generate a couple other files. One is the novoindex file (http://novocraft.com/wiki/tiki-index...ference+Genome). The others are for use with GATK (http://www.broadinstitute.org/gsa/wi...ference_genome).

3. Then, align your data using novoalign, remove duplicates with Picard, realign around indels and do base quality recalibration if desired.

4. Then, call SNPs with mpileup.

I can give more details if desired. If you want to check how mpileup is progressing, assuming you are using a command such as this:

samtools-0.1.18/samtools mpileup -AB -ugf [reference_seq.fa] [aligned_file.bam] | samtools-0.1.18/bcftools/bcftools view -bvcg -> [intermediate_file] &

Go ahead and use the following command before the first one finishes. If your reference sequence and aligned files are sorted, then looking at the end of the following file that is generated will let you know how far the variant calling has progressed.

samtools-0.1.18/bcftools/bcftools view [intermediate_file] | samtools-0.1.18/bcftools/vcfutils.pl varFilter -D99999 >[list_of_variants] &

[PratikC]

I actually haven't done this for awhile so hopefully I'm not forgetting anything.

1. Take all of your targeted regions, extend each region by X base pairs, where X is determined on your method of generating data. Put all of this in fasta format.

2. When you have your fasta file, you will want to generate a couple other files. One is the novoindex file (http://novocraft.com/wiki/tiki-index...ference+Genome). The others are for use with GATK (http://www.broadinstitute.org/gsa/wi...ference_genome).

3. Then, align your data using novoalign, remove duplicates with Picard, realign around indels and do base quality recalibration if desired.

4. Then, call SNPs with mpileup.

I can give more details if desired. If you want to check how mpileup is progressing, assuming you are using a command such as this:

samtools-0.1.18/samtools mpileup -AB -ugf [reference_seq.fa] [aligned_file.bam] | samtools-0.1.18/bcftools/bcftools view -bvcg -> [intermediate_file] &

Go ahead and use the following command before the first one finishes. If your reference sequence and aligned files are sorted, then looking at the end of the following file that is generated will let you know how far the variant calling has progressed.

samtools-0.1.18/bcftools/bcftools view [intermediate_file] | samtools-0.1.18/bcftools/vcfutils.pl varFilter -D99999 >[list_of_variants] &

Thanks Heisman, I got your point. I am doing that part and going fine with it. My question was like when we do dbSNP annotation, if our seq is mapped to ref genome, program will annotate dbSNPs on the basis of location on genome. But when we have seq reads alligned to custom reference, than how do we use dbSNP annotations?
By the way, really thans for the reply. Now I got to know how much mpileup process is completed. Can you suggest your or some other articles with similar strategy.

Thanks.

[Heisman]

Oh, that's a different question, haha.

When you create your reference sequence and make your fasta headers, I suggest making the header the genomic coordinates of the specific interval. Then, if one portion of your reference sequence is named "chr1:45673-58472" and you call a SNP at:

chr1:45673-58472 564

You can know the mutation is at chr1: (45673 + 564 -1) (be careful not to be off by one). So, even if your programming skills aren't good, you can put this in excel and convert everything to the genomic coordinates.

[husamia]

when you said you have pooled samples, I was suprised. I thought you were supposed to demultiplex your pooled samples before aligning to any reference. Otherwise you don't know which sample is which.

on another note, Is it possible you got a never ending loop in the command where you got somehow the input as your output from another process and it keeps going? just a thought like

[PratikC]

Cheers!!!

Oh, that's a different question, haha.

When you create your reference sequence and make your fasta headers, I suggest making the header the genomic coordinates of the specific interval. Then, if one portion of your reference sequence is named "chr1:45673-58472" and you call a SNP at:

chr1:45673-58472 564

You can know the mutation is at chr1: (45673 + 564 -1) (be careful not to be off by one). So, even if your programming skills aren't good, you can put this in excel and convert everything to the genomic coordinates.

Thanks Heisman, I got your point and will do as you said.
Thank you very much.
Again, can you suggest any article of yours or others showing similar kind of custom reference in NGS analysis...

[PratikC]

when you said you have pooled samples, I was suprised. I thought you were supposed to demultiplex your pooled samples before aligning to any reference. Otherwise you don't know which sample is which.

on another note, Is it possible you got a never ending loop in the command where you got somehow the input as your output from another process and it keeps going? just a thought like

Hi Husamia, yes its pooled sample and we are not demultiplexing before alignment. We just want to find out which mutation is there in population, we don't care which mutated gene coming from which patient. Well that is for just this dataset which is a pilot study. In the following main sequencing work, we will do specific patient's sequencing individually with more samples and more number of genes.

Now coming to the never ending loop, I checked that and sure there is no loop, I checked the samtools output as described by Heisman and its working properly...

[Heisman]

Here's a paper: http://genome.cshlp.org/content/20/12/1711.full

They use a positive and negative control for their pooled sequencing analysis. When we do sequencing, Phi-X is spiked into every lane, and that can be used as a negative control. A positive control is not actually necessary to run SPLINTER although it is useful for downstream analysis (similary, a negative control isn't technically necessary as you could artificially create negative control data... I doubt anybody has done that but if you are willing to miss out on the more rare variants I think you could make it work).

[aniruddha.otago]

Hi Guys,

Recently we have published a paper in Nucleic Acids research on DNA methylation data analysis- "Comparison of alignment software for genome-wide bisulphite sequence data". you might find it useful to have a look at it.. http://www.ncbi.nlm.nih.gov/pubmed/22344695

[greekkey]

I am just looking for the usage of command and saw your data. I'm wondering what's is your final result. I think you want to know the allele frequency of a given SNP in a population. If this is a common SNP (eg. 50%), probably OK, but if the freq. is not so high, in the pooled population, the signal might be masked by sequencing background noise.

Hi Husamia, yes its pooled sample and we are not demultiplexing before alignment. We just want to find out which mutation is there in population, we don't care which mutated gene coming from which patient. Well that is for just this dataset which is a pilot study. In the following main sequencing work, we will do specific patient's sequencing individually with more samples and more number of genes.

Now coming to the never ending loop, I checked that and sure there is no loop, I checked the samtools output as described by Heisman and its working properly...