Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation.
Pubmed Direct Link
Similar to another paper I blogged about, this paper describes a novel method to analyze NGS data wrapped in a biology coating (mouse genome analysis). I really like this way of describing a method since this highlights the motivation for the method's creation: to answer a question in biology.
This paper describes a method called Cufflinks, which uses TopHat to find splice junctions, which uses bowtie for mapping. Here you can see the recursive dependency tree on which this method is built. You must trust the upstream methods first.
Off to the methods section. Just like any nature paper, the methods are hidden in in the supplemental materials. I am glad to see the supplementary methods written in latex using the hyperref package. This makes it easy to move around the document quickly, plus latex is much easier to write and format than Word or other programs.
The main contribution is the transcript abundance calculation. I would strongly recommend reading the supplementary methods as this is a great Bayesian discussion of transcript assembly and abundance estimation. It also covers some nice computer science Theorems. The basic idea is to find the minimal set of transcripts that explain the fragments, then quantify each transcript. Seems simple right?
Anyhow, take a look at this paper as Cufflinks (with TopHat (with bowtie)) is a very popular pipeline for RNA-SEQ. You should probably understand what it is doing before you use it. Kudos to the authors.
Monday, May 10, 2010
Sunday, May 9, 2010
BS Seeker: precise mapping for bisulfite sequencing
BS Seeker: precise mapping for bisulfite sequencing
Pubmed Direct Link
Yet another alignment algorithm (YAAA), this time for bisulfite sequencing. This algorithm, developed in the same lab as the BS-SEQ (Jacobsen at UCLA) and called BS-SEEKER, purports to be more versatile and faster. I am of the philosophy that an aligner should be "as fast as it needs to be" to achieve the desired level of sensitivity, so the latter claim is not as interesting as the former claim. Basically, it needs to find what you are looking for first, then lets worry about speed. The abstract claims sensitivity and accuracy, so lets see if it delivers.
Surprisingly, the authors do not implement their own algorithm, but instead use bowtie. I do not recommend bowtie since it does not handle insertions or deletions, which can significantly impact mapping accuracy (see this partial discussion). Nevertheless, bowtie should be able to be replaced by other algorithms, since the novel contribution of BS-SEEKER is the separation of the four types of BS reads.
The evaluations were very lite on information, mapping reads to only human chromosome 21, and not accounting for higher sequencing error rates or indels. Furthermore, BS-SEEKER is not the best algorithm in every situation according to their own results, in some cases MAQ or RMAP performing better in terms of accuracy/sensitivity/timing. Their novel contribution is the application of bowtie, but beyond that I am not seeing the value of this method, with the exception if you use the Cokus et al. library construction. More evaluation of this method is necessary to convince others to use it.
Pubmed Direct Link
Yet another alignment algorithm (YAAA), this time for bisulfite sequencing. This algorithm, developed in the same lab as the BS-SEQ (Jacobsen at UCLA) and called BS-SEEKER, purports to be more versatile and faster. I am of the philosophy that an aligner should be "as fast as it needs to be" to achieve the desired level of sensitivity, so the latter claim is not as interesting as the former claim. Basically, it needs to find what you are looking for first, then lets worry about speed. The abstract claims sensitivity and accuracy, so lets see if it delivers.
Surprisingly, the authors do not implement their own algorithm, but instead use bowtie. I do not recommend bowtie since it does not handle insertions or deletions, which can significantly impact mapping accuracy (see this partial discussion). Nevertheless, bowtie should be able to be replaced by other algorithms, since the novel contribution of BS-SEEKER is the separation of the four types of BS reads.
The evaluations were very lite on information, mapping reads to only human chromosome 21, and not accounting for higher sequencing error rates or indels. Furthermore, BS-SEEKER is not the best algorithm in every situation according to their own results, in some cases MAQ or RMAP performing better in terms of accuracy/sensitivity/timing. Their novel contribution is the application of bowtie, but beyond that I am not seeing the value of this method, with the exception if you use the Cokus et al. library construction. More evaluation of this method is necessary to convince others to use it.
Saturday, May 8, 2010
Genome-wide mapping and assembly of structural variant breakpoints in the mouse genome
Genome-wide mapping and assembly of structural variant breakpoints in the mouse genome.
Pubmed Direct Link
This paper presents a new method to detect structural variant (SV) breakpoints called HYDRA. Since it is a Genome Research (GR) paper, the paper focuses on answering a question in biology, namely detecting SVs in the mouse genome; I will focus on the method.
SVs are detected by clustering discordant matepairs. Matepairs are pairs of reads from the same DNA fragment with some base distance (insert size) between the ends. Discordant reads are reads for which the insert size is outside the expected distribution. Hydra seeks to use multiple mappings per read, where each read may map to multiple locations, albeit with possibly different quality/probability. The method then seeks to select the correct mapping for each read collectively for all reads parsimoniously. This is to detect SVs in regions with segmental duplications, copy number changes, or other repetitive or difficult to map regions. This is an extremely important feature. A minor point is how to store multiple mappings, as a whole genome human resequencing BAM file may be hundreds of gigabytes when there is only one mapping per read. Nevertheless, HYDRA does not report what type of SV occurs, but only that there exists a SV. Further analysis needs to be performed to determine the type of SV (insert/deletion/CNV/inversion/etc.). The authors suggest to use bedtools to do this (their own software). It is yet to be seen if author software will come out that automate this whole process
Pubmed Direct Link
This paper presents a new method to detect structural variant (SV) breakpoints called HYDRA. Since it is a Genome Research (GR) paper, the paper focuses on answering a question in biology, namely detecting SVs in the mouse genome; I will focus on the method.
SVs are detected by clustering discordant matepairs. Matepairs are pairs of reads from the same DNA fragment with some base distance (insert size) between the ends. Discordant reads are reads for which the insert size is outside the expected distribution. Hydra seeks to use multiple mappings per read, where each read may map to multiple locations, albeit with possibly different quality/probability. The method then seeks to select the correct mapping for each read collectively for all reads parsimoniously. This is to detect SVs in regions with segmental duplications, copy number changes, or other repetitive or difficult to map regions. This is an extremely important feature. A minor point is how to store multiple mappings, as a whole genome human resequencing BAM file may be hundreds of gigabytes when there is only one mapping per read. Nevertheless, HYDRA does not report what type of SV occurs, but only that there exists a SV. Further analysis needs to be performed to determine the type of SV (insert/deletion/CNV/inversion/etc.). The authors suggest to use bedtools to do this (their own software). It is yet to be seen if author software will come out that automate this whole process
Wednesday, April 21, 2010
Next-Generation Instrument Google Map
Next-Generation Instrument Google Map
http://map.seqanswers.com
A map of NGS Instruments from around the world can be found my following the link above. This map is far from complete, so if you know of or have not included a sequencing facility, please add it to this map. It is amazing the number of sequencers in the field compared to only a few years ago. It is also surprising that given this global sequencing capacity, less than 20 whole genome human genomes have been published. There are an infinite number of other fun projects that do not involve human resequencing, but would it not be fun to truly have the 1000 Genomes Project completed to a high depth (more than 60X)?
I would say it is a pretty poor showing from Africa and Antartica, although BGI and BIG in China are not well reported (128 Illumina HighSeqs anyone?). The one sequencer farther from any other seems to be the one located at the "
Institute for Chemical Biology and Fundamental Medicine".
http://map.seqanswers.com
A map of NGS Instruments from around the world can be found my following the link above. This map is far from complete, so if you know of or have not included a sequencing facility, please add it to this map. It is amazing the number of sequencers in the field compared to only a few years ago. It is also surprising that given this global sequencing capacity, less than 20 whole genome human genomes have been published. There are an infinite number of other fun projects that do not involve human resequencing, but would it not be fun to truly have the 1000 Genomes Project completed to a high depth (more than 60X)?
I would say it is a pretty poor showing from Africa and Antartica, although BGI and BIG in China are not well reported (128 Illumina HighSeqs anyone?). The one sequencer farther from any other seems to be the one located at the "
Institute for Chemical Biology and Fundamental Medicine".
Thursday, April 15, 2010
High quality SNP calling using Illumina data at shallow coverage
High quality SNP calling using Illumina data at shallow coverage
Pubmed Direct Link
Here's an interesting paper on an integrated approach to SNP calling in NGS data. It is integrated in that they present both an aligner and SNP caller co-optimized for lower false positive SNP calls. It is designed for the Illumina platform and uses the four-channel (base) probability scores to better inform alignment. Where have we see that before? There are two steps, alignment then SNP calling.
The idea behind the alignment is to use the first 31 bases of the read to seed the alignment. This is based on the observation that the first 31 bases are the most reliable in any sequencing run. The only problem is that if an insertion (or deletion) occurs within these first 31 bases, then it makes it considerably more difficult to map. If an indel occurs within the first 31 bases, then the rest of the read (not the first 31) could be used to map the read. If indels are tolerated (what length indel) are tolerated in the "seed", like BWA, then no problem. Unfortunately Slider does not support indels, so you will have to re-align the data anyways if you are interested in this common source of variation. Nevertheless, their merge sort technique is a unique and novel technique among short read aligners. I'll let you read the original paper.
Their SNP caller is more interesting. It allows for priors in the form of SNP annotation databases, identifying suspect SNPs that could result from paralogs or SVs, using the per-base probability for each base in the read (so much for base calling), and the expected polymorphism rate (or homology between the reference and the current genome). I would have liked to see other tools beside MAQ being compared, as many other aligners (BFAST/BWA/SHRiMP/Novoalign etc.) are highly sensitive, while there are also a breadth of SNP callers (SAMTools/SOAP/GATK etc.). It is also unclear if the improvement demonstrated in the paper is due to the better SNP calling algorithm/model, or the aligner, or both. An easy way to test this is to substitute another aligners/SNP-caller. After filtering out low coverage regions, they also show how a SNP caller can improve concordance by using known SNPs. Be careful which version of dbSNP you use, since dbSNP 130 includes a lot of crud from NGS data projects and cancer (wont the whole genome and all alleles be present soon in dbSNP if we sequence enough cancer?).
Figure 4 is interesting since it shows that SNPs called towards the ends of reads are less trustworthy. Also, if there is high coverage at a SNP and the SNPs still occur towards the ends of reads, the concordance is extremely low (not likely anyways). Using their position in the reads seems warranted. Nevertheless, not aligning sensitive to indels will cause false-positive SNPs around where the indel should be placed, which is somewhat mediated by hypothesizing that not aligning with indels will cause false SNPs at the ends of reads. I am inclined to agree with this hypotheses from my own observations. As for the false positive reduction claim, some anecdotal evidence is suggested by Figure 6 and 7, but simulations would bolster this claim significantly. Not only does simulation provide an easy check of the claims, but also helps debug the algorithm anyways, so perform simulations when possible.
Now they need to apply this to SOLiD, add support for indels, modularize the aligner and SNP caller to allow the input/output to/from other tools, and they have an intriguing software that could be extremely useful for NGS data.
Pubmed Direct Link
Here's an interesting paper on an integrated approach to SNP calling in NGS data. It is integrated in that they present both an aligner and SNP caller co-optimized for lower false positive SNP calls. It is designed for the Illumina platform and uses the four-channel (base) probability scores to better inform alignment. Where have we see that before? There are two steps, alignment then SNP calling.
The idea behind the alignment is to use the first 31 bases of the read to seed the alignment. This is based on the observation that the first 31 bases are the most reliable in any sequencing run. The only problem is that if an insertion (or deletion) occurs within these first 31 bases, then it makes it considerably more difficult to map. If an indel occurs within the first 31 bases, then the rest of the read (not the first 31) could be used to map the read. If indels are tolerated (what length indel) are tolerated in the "seed", like BWA, then no problem. Unfortunately Slider does not support indels, so you will have to re-align the data anyways if you are interested in this common source of variation. Nevertheless, their merge sort technique is a unique and novel technique among short read aligners. I'll let you read the original paper.
Their SNP caller is more interesting. It allows for priors in the form of SNP annotation databases, identifying suspect SNPs that could result from paralogs or SVs, using the per-base probability for each base in the read (so much for base calling), and the expected polymorphism rate (or homology between the reference and the current genome). I would have liked to see other tools beside MAQ being compared, as many other aligners (BFAST/BWA/SHRiMP/Novoalign etc.) are highly sensitive, while there are also a breadth of SNP callers (SAMTools/SOAP/GATK etc.). It is also unclear if the improvement demonstrated in the paper is due to the better SNP calling algorithm/model, or the aligner, or both. An easy way to test this is to substitute another aligners/SNP-caller. After filtering out low coverage regions, they also show how a SNP caller can improve concordance by using known SNPs. Be careful which version of dbSNP you use, since dbSNP 130 includes a lot of crud from NGS data projects and cancer (wont the whole genome and all alleles be present soon in dbSNP if we sequence enough cancer?).
Figure 4 is interesting since it shows that SNPs called towards the ends of reads are less trustworthy. Also, if there is high coverage at a SNP and the SNPs still occur towards the ends of reads, the concordance is extremely low (not likely anyways). Using their position in the reads seems warranted. Nevertheless, not aligning sensitive to indels will cause false-positive SNPs around where the indel should be placed, which is somewhat mediated by hypothesizing that not aligning with indels will cause false SNPs at the ends of reads. I am inclined to agree with this hypotheses from my own observations. As for the false positive reduction claim, some anecdotal evidence is suggested by Figure 6 and 7, but simulations would bolster this claim significantly. Not only does simulation provide an easy check of the claims, but also helps debug the algorithm anyways, so perform simulations when possible.
Now they need to apply this to SOLiD, add support for indels, modularize the aligner and SNP caller to allow the input/output to/from other tools, and they have an intriguing software that could be extremely useful for NGS data.
Subscribe to:
Posts (Atom)
