今天接着看paper,突然想把以前的ChlP-Seq工作总结一下。ChlP-Seq前期或者基本的dataanalysis主要分两部分一是readsalignment,因为测序得到的read序列并不知道其在对应genome上的位置，也就是说不知道测序iT）来的read定位在genome±的什么地方，因此，首先得用alignmenttool把这些readmap到基因组上。那是不是一般blast软件都可以完成了？答案是否定的。read数目非常多，都是按照million数量级计算，并且长度短，一般为20〜3Obp左右，一般的blast软件遇到短序列，无法使用，像苦名的blasto我曾经试过，在我们实验室口己的服务器上用blat（这个可以blast短序列）mapread,消耗时间很长，最后我无法忍受停掉了，这里我不太记得我花了多少天。述有儿个问题，就是blast的吋候是否允许错配的问题。我曾经在毕业答辩的时候被问到这个问题，为什么在blast的时候耍允许错配？虽然问题很白痴，但是还是值得仔细思考。首先，我们使用的genomesequecne木身就是测序得到的，这些sequence木身就可能含有测序错谋。另外，ChIP-Seq实验小使用的样本，其sequence可能有差异，比如SNP,也就是说个体和个体直接的sequenceinformation是有差界的，并不是100%相同。还有，可能是比较重耍的一点，就是ChlP-Seq实验在sequence过程中，可能有错误。我曾经问过做ChlP-Seq实验的人（他们自己sequence序列，不是公司sequence）,ChlP-Seq实验过程中哪些因素会导致sequence错谋？其实，世界上现在对于ChlP-Seq原理并不是100%了解，尤其是ChlP-Seq实验过程中出现的各种奇怪现象的原因，人们只能在后续分析中尽量减少这些因素的影响。那么，允许儿个错配比较合适呢？目前已发表的paper±來看，都是允许2个mism珀ch。但是，没有哪一篇解释为什么是2,而不是3,4或其它。我想，可能是第一篇ChlP-Seqpaper使用的是2,于是后面的人都纷纷使用2mismatcho那如果read长度不同了？都使用2mismatch吗？这个问题值得仔细思考一下。另外，在map的过程中，只保留unique的read。为什么这样呢？因为一个read如果能map到多个位点，我们就不知道这个read信号到底是属于哪个位置？比如对于研究TF问题，我们就不知道这个TF到底是binding哪个位点。因此对于这样的read应该去掉。但是，这样去掉后，会损失很多read,我的经验是20〜30%,这个也得看具体数据。最近我也在思考，能否讣这样的read发挥它们的余热。下血说说，目前能做readalignment的比较好用的几个tool1.ELAND这个当仁不让是这方而的老大哥啦，它是Illumina公司口己开发的一个软件，速度非常快，精度也很高，这个软件我使用过，把3m订lion左右的readmap到humangenome上人概只要2小时左右，并且对内存要求很小（这个我后面会谈到）。唯一缺点就是就是和solexa测序仪捆绑销售，我想没人会发疯到，为了用ELAND而去买台solexa测序仪吧？这个软件低版本最长只支持到32bp的read,新版本ELAND解决了这个问题。2.SOAP这个是ELAND很好的替代品，而且是中国人写的，其中还有一个小孩，是北大的，很NB。平时不管我问他什么问题，他都或多或少懂一些，并H.能跟我讲的很清楚。这个软件速度也比较快，但是比ELAND慢，前面同样的数据，SOAP人概需要1天对一点时间才能完成。这个软件很要命的一点就是，对内存要求很高。因为它是把基因组信息读到内存中建索引，大概是genomefile大小的4倍！如果要做human的map,至少需耍32G内存，一般小型服务器是无法满足的。这个软件有个参数很有意思，大致就是，考虑到测序过程中，对于一个read测得序列越长那么出错的概率就会越人，于是在map过程中，就依次截掉一个3'碱基，然后再map,直到序列太短。当然这里面具体的参数可以设置。此外，SOAP在map过程中允许gap,这是ELAND做不到的。如我前而所谈到的，我很想弄清楚，这些sequencingerror是什么，但是冃前述不知道。具体可以看Li,R.,Li,Y.,Kristiansen,K・Wang,J.SOAP:shortoligonucleotidealignmentprogram.Bioinformatics24,713・714(2008).3.Maq这个也是中国人写的，这个人似乎现在在英国。我不太喜欢用这个软件，因为它用一个mapquality参数去衡量map的结果，如果我想要unique的read,这个软件就比较烦了。当...