31

从WGS测序得到的VCF文件里面提取位于EXON区域的

首先要下载并且得到人类基因组的外显子坐标记录文件

这里我用的参考基因组版本仍然是hg19,所以去CCDS数据库里面下载对应版本,并且格式化成BED文件。

wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/archive/Hs37.3/CCDS.20110907.txt
cat CCDS.20110907.txt |perl -alne '{/[(.*?)]/;next unless $1;$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$" foreach split/,/,$exons;}' >hg19exon.bed

Continue reading

19

测序深度和GC含量的问题!【直播】我的基因组 47

在前面我们提到了用ChIP-seq的分析方法可视化了一下我的WGS数据,结果我们的测序深度分布居然是跟基因组的genomic feature相关的~~~
比如在TSS附近,就很明显看到了一个测序深度峰值,那么前面我们并没有给出直接的解答,而且简单的提了一下这是二代测序的特点,GC含量片段偏好性!
作为一个合格的生物信息学工程师,我当然要把这个理论用自己的代码和数据来亲身实践一遍~
我首先把全基因组的bam文件用mpileup模式输出,根据1000bp的窗口滑动来统计每个窗口的测到的碱基数,GC碱基数,测序总深度!
代码比较复杂,一般人可能理解不来的!

Continue reading

十二 19

【直播】我的基因组22:用IGV查看具体某个位点是否变异

下载IGV和导入文件的方法我就不多说了,可以直接在windows平台下使用,就跟你操作QQ一样,自己摸索就好了!

著名芬兰运动员Eero Mäntyranta,他拿过七枚奥运奖牌。他的血红细胞远超正常人水平,甚至一度被奥组委误以为服用了禁药。后来经过研究发现,他的EPOR基因上的一个位点rs121918116,发生了一个G>A的变异,使得他的血氧含量达到了普通人的150%,所以他耐力惊人。

在snpPedia里面可以查看这个位点的信息:http://snpedia.com/index.php/Rs121918116 Continue reading

十二 19

【直播】我的基因组(20):覆盖度详细探究

前面我们在第8讲提到了公司给我的一个报告的统计表格,有人反映不会做。

本来应该只需要给我6亿条reads的(PE150测序,人30X),但是足足给了我8.9亿条!(但事实上很多paper发表的基因组高于60X的也不少)

表格里面提到了好几个概念,比如duplicate的reads,一般说的PCR造成的duplicate,在找变异的时候需要去除掉。然后是那些比对到了不同染色体的reads pair,虽然只有2.29% ,也是需要重点分析的。(前面我也讲了如何提取以及分别分析它们!) Continue reading

十二 19

【直播】我的基因组(18):初步分析PCR duplication的情况

我博客里面有详细讲读原文查解去除PCR duplication的reads的原理和方法,还比较了samtools和picard这两个软件的区别,请点击阅看(仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的),或者复制链接(http://www.bio-info-trainee.com/2003.html)到浏览器查看。 Continue reading

十二 19

【直播】我的基因组(17):初步分析一下multiple mapping 的情况

分析multiple mapping的情况,还是先用公司已经提供的bam文件来吧,后面我会用自己比对得到的bam文件再重新分析一次。

公司给我的bam数据是把有效测序数据通过 BWA(Li H et al.)比对到参考基因组 (b37),这是目前使用率很高的一个软件。在比对过程中,如果一个或一对 read(s) 在基因上可以有多个比对位置,BWA 的处理策略是从中选择一个最好的,如果有两个或以上最好的比对位置,则从中随机选择一个。这种多重比对(multiple hits)的处理对 SNP、INDEL 以及 CNV 等的检测有重要影响。 Continue reading

十二 19

【直播】我的基因组(16):提取左右端测序数据比对到不同染色体的PE reads

这类情况仅仅针对于双端测序数据,因为根据实验原理来看,对一个DNA片段,会把它的左右两端分别测序,但是测序仪器的测序长度有限,对本次实验来说,打断的DNA片段长度在350个碱基左右(这个长度只是一个分布,并不是真实值),理论来说测序是左右各150,加起来也就300,也就是说DNA片段中间还有50个碱基是测不到的(当然,实际上是有可能测通的)。而对这个配对的reads来说,来自于同一个DNA片段,所以理论上它们应该比对到同一条染色体的。也还是基于对sam格式的文件的理解,前面我们提到了sam文件的第3,7列指明了该reads比对到哪条染色体,以及该reads的配对reads比对到了哪条染色体(如果比对到同一条染色体,那么第7列是=符号)。所以我们只需要写脚本来提取即可! Continue reading

十二 19

【直播】我的基因组(15):提取未比对的测序数据

之前我们说了比对上的数据,那么会有人想到有没有没有比对上的数据呢?

既然是从我的血液里面提取到的DNA进行测序的,那么理论上测序仪出来的所有测序reads都应该是我的,也应该都可以比对到人类的参考基因组,但是实际过程中的确存在着未必对上的数据。

现有人类基因组毕竟只是个参考,也许我有某些独特的DNA序列呢?而且也不一定测序数据就都是人类的,也许我血液里面会有那么些微不纯粹呢?也有可能是某些片段变异的太多了,超过了常见的比对软件的承受能力,可以试用SHRIMP这个软件来提取一下。当然,这不是本讲的重点。 Continue reading

十二 19

【直播】我的基因组(14):bam文件给按照染色体给分割成小文件

昨天,我们了解了一下SAM格式的比对结果,不知道大家理解的怎么样。但是全基因组测序数据实在是太大了,即使比对后把sam文件压缩成二进制的bam文件也还有55G(如何压缩转换可查看直播十二),如果完整的导入IGV查看会略微考验计算机配置。

Continue reading

十二 19

【直播】我的基因组(13):了解sam格式比对结果

十一讲中将我们主要讲了如何将下机数据比对到参考基因组中。但是很多人对比对结果却是一头雾水。那我们现在来了解一下Sam格式的比对结果吧!

比对工具到现在已经多如牛毛了,见列表: https://en.wikipedia.org/wiki/List_of_sequence_alignment_software 。但是能被大多数人熟知的,就是bowtie和bwa(我们在十一讲中用的才是bwa),它们把测序数据比对到参考基因组之后,都会生成一个sam格式的文件。随后的大部分分析都是基于sam格式进行的分析,虽然Jimmy多次强调这些基础知识的重要性需要大家私下自学。但是由于这个sam文件实在是太重要了!!!所以,不得不亲自抽出一讲来说说它,后面也会基于此写十多篇文章: Continue reading

十二 09

【直播】我的基因组(十二):先粗略看看几个基因吧

昨天我们说到,测序得到的fastq文件map到基因组之后,我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format,而BAM就是SAM的二进制文件。通常sam文件太大,我们会生成bam文件来节省空间。sam文件和bam文件的转换用samtools这个软件就可以完成。 Continue reading

十二 09

【直播】我的基因组(十):测序数据质量控制

质控之前我们在直播八的时候分析过,公司也给了我质控后的的数据,但是毕竟是别人做的,我们做为一个数据分析师,自己动手来验证一下公司给出的报告也是再好不过的了。大家可以跟着我先将下载数据进行一下质控。 Continue reading

十二 09

【直播】我的基因组(十一):测序数据的比对

上一次直播中,我们对拿到手的测序数据进行了质控,测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了,这里我们先用参考基因组hg19,大家可以参照【直播】我的基因组(五):测试数据及参考基因组的准备来下载参考基因组hg19,我这里选择的是UCSC提供的hg19。然后安装bwa软件进行比对,可以参考【直播】我的基因组(四):计算资源的准备来安装,以及对hg19建立索引。 Continue reading

十二 09

【直播】我的基因组(九):拿到数据后要做的事情

时隔好几个月,因为各种各样的原因数据终于拿到了自己的手上,真是不容易啊!

拿到数据后,第一件要做的事情就是检查数据传输的完整性,然后备份!我拿到的数据如下:

可以看到,公司给了我测序仪的下机数据(raw data)和他们质控后的clean data,这个过程减少了6G的数据量,对应着约90亿bp的碱基,相当于减少了3个人的全基因组数据。具体推算公式见前面的系列直播贴!

首先我把数据拷贝到了我上上周买的2T移动硬盘里面,再拷贝到我工作电脑一份,服务器一份,私人电脑一份,另外一个移动硬盘一份。然后删除了公司寄给我的硬盘里面的数据,再把硬盘寄回给公司,然后监督他们删除我所有的数据。(做这么多就是为了保护隐私,当然这个大前提是我已经确定数据没有问题了。)

检查数据传输的完整性就是md5校验,看看数据在拷贝过程中有没有意外的损坏(这个在之前下载数据的时候我也说过)!一般传输数据之前,会用md5命令来生成各个文件的md5值,就是下面的MD5.txt文件里面的内容,然后传输数据之后,需要自行用md5sum -c MD5.txt 来校验文件里面记录的文件的完整性,如果显示都是OK,说明文件拷贝传输过程是没有问题的!但这个过程会耗费大量的磁盘读写,磁盘读写能力是有限的,所以开多个进程并不能加快这一过程。

然后我把公司处理好的bam文件上传到服务器做下游分析,我用的winscp软件把文件传到服务器上的!

从明天起,我们就开始正式对基因组进行分析啦!欢迎围观!

请扫描以下二维码关注我们,获取直播系列的所有帖子!

1

十一 23

【直播】我的基因组(八):原始测序数据质量报告

由于我是分期付款,所以我先拿到了我的测序数据的质控结果和比对情况分析报告,需要补齐全款后才能拿到原始测序数据!(中间还出了个小意外,打款的时候不小心多打了30块钱!(⊙o⊙)…不过多打的30块钱想拿回来估计不太可能了,需要填写书面申请表格并且自费快递到公司,这边跨境快递费都不止这个数了) Continue reading

十一 23

【直播】我的基因组(七):从整体理解全基因组测序数据的变异位点

首先记住一个很重要的知识点,变异是相对的!

简单说一下什么是找变异,变异跟突变有什么区别呢?举个栗子:有国际组织规定了人类的参考基因组(如UCSC,ENSEMBL,NCBI等,前面帖子都有讲),就是 AAAAA(这里简化一下,就5个碱基,其实人类基因组多达30亿个)  。现在通过给自己测序得知,我与之对应的是AGCAA,那么我相比国际基因组来说,就是2个变异位点,位于基因组的坐标2和3,但是它们还不能说就是突变。 Continue reading

十一 23

【直播】我的基因组(六):变异位点注释数据库的准备

通常一个人的全基因组测序数据可以挖掘到四百万个SNVs(跟参考基因组不一样的单碱基位点),还有五十万的indels(insertions or deletions),但是得到的数据通常是以vcf文件格式给出的(自行搜索什么是vcf格式),比如下面:

很明显,正常人是看不懂这些变异位点有啥子一样的,只知道第20条染色体的1230237坐标上面本来是一个T碱基的,但是突变成了G,那么我们必然还想知道,这个位点是在某个基因上面吗?如果是,在基因的外显子还是内含子?它的突变有没有改变该基因的功能呢?有没有影响它的转录和翻译呢?还有世界上有没有其他正常人也是这个位点变异呢?如果有,是哪些人种呢?有没有癌症病人也发现了这个变异呢?如果有,是什么癌症呢?所以我们必须下载一系列的变异位点注释数据库,来全方位的解释我们自己找到那四百万个SNVs和五十万的indels。下面我们一起进行数据库准备。

Continue reading

十一 07

【直播】我的基因组(五):测试数据及参考基因组的准备

我的全基因组数据还没拿到,而且还会推迟,简单说(tu)明(cao)一下原因(还好当初为了避免广告嫌疑一直没说是哪个公司负责测序,反正用的是illumina的hiseqX10这个测序啦,所以可以尽情的吐槽)。

心烦意乱的吐槽线 Continue reading

26

【直播】我的基因组(四):计算资源的准备

大家久等了,Jimmy的测序数据还没有拿到手!但是,工欲善其事必先利其器!所以jimmy在等待自己基因组的这段时间里,准备好了自己计算资源!鉴于会有不少同志们会跟着直播来自己动手分析一个公共的全基因组重测序数据(参考数据下载地址详见直播一),我在这里和大家分享一下计算机资源及软件的基本要求。 Continue reading