序列比对与NGS

序列比对

NGS

43 thoughts on “序列比对与NGS

  1. 博主好,有个测序的问题想请你帮忙解答一下。
    我刚做了一个植物不同组织不同处理条件下的RNA-SEQ,得到了很多数据。其中有一项是统计已得到转录本的RPKM值及分布区间,我不解的是,我的数据的RPKM大多在0~1区间,超过55%;然后1~3,3~15,15~60区间的都在10%-15%左右,RPKM>60的只占5%。我想知道这样的分布是否不正常,原因可能是什么?
    另外,还有个问题,我做的这个测序是有参RNA-SEQ,但我想知道测序深度和覆盖度,这个怎么计算?测序公司说有参转录组测序一般不计算深度和覆盖度,那我怎么判断测序结果覆盖了多少基因组?
    谢谢,希望得到你的指点。

    1. 测序结果应该算正常,rpkm值应该符合power law分布原则,你可以把你的结果画一个图,一看就清楚了。但是不知道测序深度,所以无法全面判断。一般而言,有参考基因组的计算测序深度并不困难。所谓的有参不考虑测序深度显然是不精确的,它只相较于无参的而言,不需要那么高的覆盖度。你可以查看博文illumina NGS常见问题来了解一些基础知识。除了覆盖度,您还需要了解的是你的fastq文件的质量,使用fastqc就可以轻松得到。之后需要了解的是mapping的reads数,可以参考如何统计BAM文件中的reads数。祝你好运。

  2. 谢谢博主的耐心解释。我又跟测序公司沟通了表明想知道测序深度和覆盖度问题,并给他们发了一个链接http://emuch.net/html/201110/3761305.html,他们按这个链接给我算了一组数,说是覆盖度,
    CKL
    9.83421120798756
    CKR
    11.0388165124084

    他们说有参RNA-SEQ没必要纠结这两个概念。
    想进一步请教博主,有没有办法自己测算这两个指标。另外,怎么通过测序数值计算我材料的转录组大小。
    再次感谢博主帮忙。

  3. CKL和CKR分别是我两组样本的编号。
    博主能否将电邮或即时联系方式告知我,我很急切的想把一些问题解决掉以便后面的数据整理和筛选工作。我的电邮是tea_gl@126.com。谢谢!
    另外,恳请博主先帮忙解释一下上次的问题:
    有没有办法自己测算我这个RNA-SEQ的深度和覆盖度这两个指标。另外,怎么通过测序数值计算我材料的转录组大小。
    再次感谢!

  4. 你好,cufflinks可以对以sam格式储存比对结果进行组装,你的博客也提及到cufflinks可以组装tophat的输出结果,通过cufflinks说明,我们可以知道它支持别的比对软件的输出的结果,比如bwa。我目前不清楚的是不同软件输出的比对结果不是一致,尽管都是以sam格式储存,但是各个软件输出的结果在sam里面却是不一致的,比如,一个read能比对到多个位置上,在tophat输出的sam里包含着每个比对位置的信息,并且有NH:i:5类型的tag描述比对次数的,有的软件对于1个read比对到多处,只随机输出1次,另外还有的比对软件支持clip alignment,这个软件输出的sam里的‘CIGAR’对应列会有‘S’字符,据我所知tophat中的’CIGAR’是没有’S’字符的。我目前不太清楚cufflinks是否能正确处理各种软件输出与tophat不同编码比对信息形式sam文件吗,迫切的想了解清楚,因为这涉及到cufflinks能否正确对基因和转录本进行准确定量的问题?不知道你知道这方面的信息,希望不吝赐教,谢谢!

    1. 不好意思,我并不了解这方面的准确信息。建议您先阅读其主页上的文献,上面可能有您需要的答案。

  5. 博主,你好。从你的博客里学了不少东西。我从ENCODE上下载到的通过cufflinks处理后的数据(.gtf/.gff)文件,里面对应的gene id是cuff.xxxx, 我想将其map到refseq或者ensemble上,该如何映射呢?

    1. 如果你会使用R/Bioconductor的话,你可以使用由我维护的ChIPpeakAnno来完成这项任务。如果不会,你可以bedtools来完成注释。

  6. 博主,您好!
    请教一下我在使用cuffdiff的时候splicing.diff里面的结果都是不显著,结果如下
    TSS100009 XLOC_062179 – chr15:96273598-96369222 s484 s485 NOTEST 0 0 0 0 1 1 no
    TSS100010 XLOC_062180 – chr15:96273598-96369222 s484 s485 NOTEST 0 0 0 0 1 1 no
    而且value值都是0 怎么回事的呢

    1. value为0只是表示rpkm值非常小,在输出时用0表示。你需要先检查一下你的mapped的reads数,查看一下是否是测序深度不足。

  7. 博主,您好 !后来我修改了一下参数 添加参数–min-js-for-test,但是博主 我想问一下gene_exp.diff中很多差异表达的基因都是
    XLOC_001863 XLOC_001863 ENSG00000187147 chr1:44886177-44886640 s4 s5 OK 36.3725 0 -inf -nan 0.00035 0.0447368 yes
    XLOC_001864 XLOC_001864 ENSG00000187147 chr1:44887789-44888404 s4 s5 OK 26.3779 0 -inf -nan 0.0002 0.036662 yes
    总有一个表达为0,这是怎么回事的呢
    在splicing.diff中

    TSS10753 XLOC_006009 ENSG00000198771 chr1:167604197-167604377 s4 s5 OK 0 0 0.832555 0 5.00E-05 0.000578641 yes
    TSS107588 XLOC_063162 ENSG00000140650 chr16:8907368-8907537 s4 s5 OK 0 0 0.832555 0 5.00E-05 0.000578641 yes
    TSS108242 XLOC_063575 ENSG00000005189 chr16:20849735-20850656 s4 s5 OK 0 0 0.832555 0 5.00E-05 0.000578641 yes

    很多test_stat为0 而且sqrt(JS)很多很相似,这是什么原因的呢?

    1. 你是指的–min-reps-for-js-test参数吗?你有几个样本呢?这个参数应该设置为样本数的一半或者以上。出现一边是0的情况很正常。但是这些基因都不应该是主要考虑的阳性。还有就是对于单exon的转录本,也需要十分小心。

  8. 请问,利用cuffdiff进行对照组和病例组的差异表达分析时,是以什么为参考基因组啊?是以他们组装成的一套转录本呢还是以对照组的转录本或者病例组的转录本?

    1. 如果需要发现新的转录本的话,需要使用cuffmerge把所有的转录本合并起来一起使用。其它情况,可以使用已知的转录本。

  9. 博主,您好! 我主要有三个比较样本 但是均无replicate所以得设置这个参数。
    1、那么针对单exon的转录本,是不是在进行cuffdiff的时候就必须事先过滤的呢?
    2、还有博主我想问一下在做基因表达和可变剪接的时候适合在一个流程做得么?
    3、如果想识别新的isoform gene的时候,是不是在tophat比对的时候不加入注释文件会比较好的呢?
    4、cuffmerge cuffcompare都可以生成cuffdiff需要的注释文件,而且对最终识别结果数目还是有一定的影响,有什么选择的建议的呢?

    1. 1对于单exon转录本,不需要事先过滤。这不影响最后的结果。而且,这些也是有效的reads,对于估计样品总量是有帮助的。2.对于基因表达和可变剪接,我还是推荐tophat2->cufflinks->cuffmerge->cuffdiff这个流程。3.如果想识别新的isoform gene,最好还是要加入已知的注释,否则可能假阳性比例会很高。当然,你也可以都做一遍,然后比较。4.使用cuffmerge。我认为你可以仔细阅读一下cuffcompare是做什么用的。

  10. 博主你好~
    请教一个在GTF文件里标记start/stop codon的问题
    我想计算基因组中stop codon到last intron的距离分布,但是genome的gtf文件中未标出stop codon的位置。
    我现在有的是所有基因的Coding sequence,可否基于这些序列,来在gtf上标记出start,stop codon?

    感谢。

    1. 我想你可以试着使用R/Bioconductor中的GenomicFeatures来调出所有的UTRs,或者CDS。这样你就知道有哪些start, stop condon了。我觉得使用coding sequence也是可以的。但是对于一个基因有多个start, stop condon position的情况,你需要多加注意。

  11. 你好,看过你的博客之后, 对我很有帮助。 我接触RNA seq data 也有半年了, 但接触的都是mRNA的分析。 我最近收到一份small rna的数据, 是由illumina产生的, 并且有producer提供的adapter 序列。 但是我对数据使用cutadapt去除adapter之后, 大部分的reads长度峰值都在30nt, 但是理想中的峰值应该在24nt左右, 我尝试过其他的adapter序列, 并通过fastqc提供的污染序列进行trim, 峰值总是在30nt。 我看过你的关于barcode的文章, 请问我需要也对数据去除barcode么? 该怎么去除? 能同时和adapter一起去除么?还是有先后顺序? 万分感谢!

    1. 这与你的实验有关。如果你的barcode是在3’adapter中间的,那估计不需要去去除。如果你的barcode是整合在测序结果里的,就是5’端最开始的几个,则需要去除。去不去除barcode取决于你有没有加barcode和怎么加的barcode。至于为什么是30nt,似乎看上去好象是正好多了6个碱基,但是具体还是需要看实验。barcode的作用有很多,但主要还是用来区分混合在一起的多样品的来源,它可以和adapter同时去除,但是需要你设计好程序。

      1. 嗯,太感谢了,理解了。 那请问如果实验里加入了barcode, barcode会出现在fastq文件里所有的reads里么? 还是部分reads会有barcode? 还有我专门把峰值30nt的reads都输出到一个文件里, 发现~94%的reads 都是一模一样的序列, 这种情况正常么? 毕竟我是按照数据产生者提供的adapter进行trim的, 但是trim完之后的结果却没那么满意。 谢谢!

        1. 94%是同一序列有点不太正常。你先查一下这段序列是什么。应该说不可能每个reads里都有barcode,但是如果有的话,应该是大部分都有。

          1. 比如对于一个sample, 在5053220的30 nt的 reads中有4777977个序列是“TGCTTGGACTACATATGGTTGAGGGTTGTA”。 这reads里面也不包含对应的barcode. 所以我一直在处理data上面花了很多时间。 谢谢你的帮助呢, 收获很多! 很赞的博客, 会一直支持的。

          2. 这段序列好象在很多果蝇的miRNAseq的结果中都出现过。你先把它记录下来吧。我blat了一下,发现它在果蝇的基因组中出现的100%mapping的结果有很多。可能这是一个正常的现象。但我没有做过果蝇的miRNAseq分析,所以不能提供更多的帮助了。

  12. 博主你好,不知道你对bwa这个软件了解多少,如果可以的话能不能讲讲这个软件,特别是软件三个算法中需要设置的各个参数所代表的意义

  13. 博主,您好,我想统计一下经过bwa或者tophat map之后,统计reads在基因组相应区域的个数,比如在exon区有多少reads。 我已获得bam格式文件,请问有什么办法吗?

  14. 请问博主,当使用gtf为denovo组装得到的结果,基因与转录本一一对应,为什么使用cuffdiff做差异分析的时候得到的差异基因个数和差异转录本个数会不一致?

    1. 您为什么会认为它们应该一致呢?您的意思是每一个转录本都被视为一个基因?除非这些基因或者转录本完全没有重叠,否则请参见:https://www.biostars.org/p/13525/

  15. 博主你好,我想问问我们做靶外显子测序,参考基因组选择全基因组呢,还是蛋白编码转录本序列,还是全基因组过滤掉靶基因意外的染色体?GENCODE和UCSC选择哪一个呢?还是都可以呢?谢谢

  16. 老师您好,感谢您百忙之中看我的提问。
    现在遇到一个RNA_seq的解析问题。想咨询您以下:
    我想解析GSE63055这个数据,里面是paired end sample,但是,每一个sample的数据都有4个experiment,一共4组8个数据,平时见到的paired end的sample都是1个experiment,1组2个数据,使用tophat的时候就一个命令里面直接计算了。但是这回是8个数据,使用tophat的时候应该如何计算?最后出来的FPKM数据肯定是一个sample对应一个FPKM数值,这样4个experiment的话,一个sample就会出现4个FPKM数值,那是取这个4组数据得到的FPKM值的平均数作为这个sample的FPKM值么?还是在一个tophat的命令中,把这8个数据都放进去一起算?不知道我表达的清楚不清楚,还望老师百忙之中可以解答。

    1. 对你的问题,我的理解是两组,每组里有四次重复?如果是这样的话,tophat mapping之后使用cufflinks就可以了。cufflinks可以一次把多个重复的数据都传进去.

      1. 谢谢老师这么快速的回答,那个数据的样子是这样的,最终结果是1个sample对应一个FPKM值,正常的话,如果1个paired end sample里面是1个experiment,会有2个fastq数据,tophat之后得到一个bam file,cufflinks以后应该是得到一个FPKM值,应该是这样吧?而现在这个是1个sample里面,有4个重复,每个重复有2个fastq数据,这样tophat以后应该是得到4个bam file,然后将这四个bam file一起通过cufflinks计算以后得到一个FPKM值么?不好意思,问题有点长

  17. 博主您好!从您的博客学到超级多的,建议博主也开个LncRNA分析的贴啊,谢谢啦!

发表评论

电子邮件地址不会被公开。 必填项已用*标注