十二 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这个软件就可以完成。

samtools view -h abc.bam > abc.sam
samtools view -b -S abc.sam > abc.bam

我们已经拿到了bam文件,我这里就先用公司给我的bam文件吧,根据我的帖子:仅仅对感兴趣的基因call variation ,可以先了解几个比较有趣的基因的变异情况。我自己呢,对以下几个位点和基因比较感兴趣,就用他们来讲一下今天的内容吧!

1.STAT4上的rs7574865和HLA-DQ的 rs9275319是中国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因

2.V1aR基因是雄性标志性出轨基因。

3.GLI3和PAX1基因控制鼻孔的大小,而RUNX2基因控制鼻梁的宽度。DCHS2基因调控鼻子的突起程度,即决定鼻尖是否朝上和鼻尖的角度,或者说它决定了你的鼻子是否迷人挺拔。

4.肥胖有关的基因FTO(Fat Mass and Obesity Associated),最近发现了调控肥胖(主要是脂肪燃烧)的基因是IRX3 和IRX5。大约100个基因位点与BMI(身体质量指数)相关,600个基因位点与身高相关,160个基因位点与肥胖特征如腰臀比相关。6个新基因位点,这些位点位于LEMD2、CD47、GANAB、RPS6KA5/C14orf159、ANP32和ARL15基因内或周围。

那,我们就先关注这几个基因吧(不要问我为什么(-_-メ) )。

首先找到这些基因的坐标,看到如下:

其中V1aR基因这个雄性标志性出轨基因,在标准的基因命名系统里面其实是AVPR1A:http://www.genecards.org/cgi-bin/carddisp.pl?gene=AVPR1A ,这里面涉及到HUGO symbol的概念,这个genecard数据库也非常赞,基因相关信息都可以在这里面查找的。

有了这些坐标信息,我们就进入我们的基因组工作目录:

cd data/project/myGenome/

然后把坐标文件做好

因为公司给我的bam文件里面,用的参考基因组是GRCh37而不是hg19(两者区别在于chr是否标记),我们还是需要下载;

cd ~/reference

mkdir -p  genome/human_g1k_v37  && cd genome/human_g1k_v37

# http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/

nohup wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz  &

gunzip human_g1k_v37.fasta.gz

wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai

wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/README.human_g1k_v37.fasta.txt

然后回到基因组工作目录,保证bam文件在上图中bamFiles那个目录,然后用下面这个脚本,批量提取我们感兴趣的基因的变异情况:

cat key_gene.list |while read id;

do

chr=$(echo $id |cut -d" " -f 1|sed 's/chr//' )

start=$(echo $id |cut -d" " -f 2 )

end=$(echo $id |cut -d" " -f 3 )

gene=$(echo $id |cut -d" " -f 4 )

echo $chr:$start-$end  $gene

samtools mpileup -r  $chr:$start-$end   -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta bamFiles/P_jmzeng.final.bam  | \

bcftools call -vmO z -o $gene.vcf.gz

done

等三分钟就好了,结果如下:

前面我们说到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因,那么我们很容易去dbSNP数据库或者我最近强烈推荐 的snpedia数据库(吐血推荐snpedia数据库,非常丰富的snp信息记录)里面找到它的坐标。

6 32666295 :Rs9275319--HLA-DQ

2 191964633 :Rs7574865--STAT4

然后我检查了我刚才call到的variation文件,

zcat STAT4.vcf.gz |grep -w 191964633 显示为空。

zcat HLA-DQ* |grep 32666295  也是空。

哈哈,我完美的错过了这两个易感位点!!!!谢天谢地!!!

其余的我就不讲了,毕竟会涉及到隐私,我就讲这个方法吧!

文:Jimmy、吃瓜群众

图文编辑:吃瓜群众

十二 09

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

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

因为此直播系列走得是半科普半技术路线,所以我这里show一个最常用也是最简单的测序质量控制软件,大名鼎鼎的fastqc软件,它是一个java软件,功能很单一,就是对你的测序数据生成一个网页版的可视化检测报告而已。这个软件的安装可以查看之前的直播贴(【直播】我的基因组(八):原始测序数据质量报告)。它在在linux或者windows平台都可以使用。直接下载这个压缩包: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip

我比较喜欢把这些软件放在biosoft文件夹下面(个人windows电脑),这个软件安装后会自带一些数据,大家感兴趣可以查看一下。

由于fastqc是免安装软件,直接解压后就可以直接使用。解压打开里面后缀是 .bat (相对于windows平台的批处理程序)的文件就打开fastqc啦,然后导入数据开始分析即可,静候一两个小时。

如果你用的是linux服务器,可以直接用unzip解压fastqc的zip压缩文件。里面有个fastqc的文件,就是fastqc的程序了。我们可以用fastqc  -o output dir [-(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN的命令让它进行质量控制。-o是用来指定输出文件的目录,注意是这里是不能自动新建目录的。输出的结果是.zip文件,默认自动解压缩,-noextract则不解压缩。-f用来强制指定输入文件格式,默认会自动检测。-c用来指定一个contaminant文件,fastqc会把overrepresented sequences往这个contaminant文件里搜索。后面加上你要质控的序列的文件名就可以了。

把所有的fastq.gz文件用fastqc软件处理得到的测序质量检测报告是一个html文件加上一个文件夹,如果没有解压缩需要用命令ls *zip|while read id;do unzip $id;done,把所有压缩包批量解压开。可以看到对每个测序数据它都进行了十几项统计结果和可视化的图片,对该款软件的结果感兴趣的可以下载(http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip) 文件,对原始数据处理前后的fastqc报告的区别显而易见。

然后批量抓取里面的%GC,Total sequences等信息,来跟之前公司给我的报告做比较,看看公司给我的结果有什么出入!

我以前写过帖子关于如何得到fastqc的统计表格:写脚本对fastqc的结果进行统计咯!


也就是说把多个qc的结果通过脚本整理在一起,方便查看。我们的统计结果如下:

当然一般不会有什么差别的,而且fastqc跑出来的结果都是合格的,公司对raw data得到clean的步骤仅仅是过滤掉不合格的reads,全部丢弃,而不是截断,豪气!!!

因为illumina的X10机器跑出来的数据一般都非常不错,我就没有在这里面下太多功夫,只是走个流程看一下测序质量,的确非常好,大家如果遇到质量比较差的数据,可以去我博客里面寻找各种解决方案。当然,质量控制不只是看序列的质量,还有很多小技巧,我会在后面的帖子里面专项讲解,比如我的数据是5条lane的数据合并起来的,那么lane的上样品是一定正确吗,那些没有比对上的reads是什么之类的相关问题。

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

1

十二 09

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

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

我们首先简单讲一下为什么要进行比对以及比对过程是怎样的?

可以看到我们到手的测序数据格式是fastq,每条reads都是150个碱基组成,如果只看这fastq,我们没办法得知它的意义,参考基因组那么大(人类约30亿个碱基),这个reads在我们基因组的哪里呢?

简单解释一下,假设人类基因组是123456789,如果我们的测序得到的reads是123,那么我们很明显知道这条reads在基因组的第一个位置,跨越了3个长度,如果我们的reads是567,那么我们也可以看到它在基因组的第5个位置。如果我们看到了一个reads是235567,同样我们也很容易看到它在基因组第2位置,但是在第3个位置参考是4,它却是5,这里可能是测序错误,也可能是这个reads真的变异了!

但是我们的参考基因组远远没有那么简单,30亿个碱基的庞大数目,测序的一条reads也有150个碱基,仅仅用肉眼观察基本不可能判断出它到底在哪里。但并非一定观察不到,如果你有多的不可计的时间及精力的话,手工比对穷极一生来搞定一条reads的比对就很不容易了(当然肯定不会有人这么傻,这里只是说数据量真的很大而已)。然而在我们手上可是有8.9亿条reads,所以我们需要借助计算机来进行比对,现在比较流行的基因组比对工具是bwa和bowtie,它俩的算法不一样,但是我们不需要了解那么具体,只需要知道它可以把我们的fastq测序文件通过与参考基因组的比对生成sam格式(自行搜索了解该格式)的比对结果文件(如下),从sam文件中,我们可以看到每条reads在参考基因组的位置,这条reads是在哪一条染色体,又是在这条染色体的哪个位置就可以一目了然。

对于比对的结果,我们可以用IGV可视化查看,还可以手动查看每个基因的比对情况:

下面我简单讲一下代码

1,下载hg19基因组

cd ~/reference

mkdir -p genome/hg19  && cd genome/hg19

nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &

tar zvfx chromFa.tar.gz

cat *.fa > hg19.fa

rm chr*.fa

首先要理解linux基础命令,在我们的服务器上面新建好目录,找到hg19的下载链接,用linux自带的wget下载,因为文件太大,所以我们用nohup放在后台下载。下载后是压缩文件 chromFa.tar.gz,在linux里面需要用tar zvfx 来解压tar.gz文件即可。解压开后是一个个文件,需要用cat合并!最终效果如下:

2,安装bwa软件

## Download and install BWA

cd ~/biosoft

mkdir bwa &&  cd bwa

#http://sourceforge.net/projects/bio-bwa/files/

wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2

tar xvfj bwa-0.7.15.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files

cd bwa-0.7.15

make

~/biosoft/bwa/bwa-0.7.15/bwa

我所有的软件都安装在自己的home目录下面的biosoft文件夹。同样,也是找的bwa的下载地址,然后解压,然后直接make即可。很多人的服务器会报错zlib.h缺少的问题,看我以前的教程:http://www.bio-info-trainee.com/518.html ,缺少什么你就安装什么,但是缺少的东西需要安装到系统环境变量,但是我的bwa是直接安装到自己的目录,所以我用全路径在调用该软件。如果你的这个命令~/biosoft/bwa/bwa-0.7.15/bwa  能够显示下面的help文档,说明你已经安装成功啦~

3,对hg19参考基因组用bwa构建索引

cd ~/reference

mkdir -p index/bwa && cd index/bwa

nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index   -a bwtsw   -p ~/reference/index/bwa/hg19  ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1   &

代码很简单,就是新建好一个文件夹来存放我们的参考基因组的索引,我这里选择的是我的home目录下面的reference/index/bwa/ 文件夹,可以看到如下内容:

我还是用了nohup把这个命令挂在后台,防止掉线,因为要运行2个小时左右,我加上time命令可以看到运行时间,我用了bwa的index模式来索引参考基因在,具体bwa用法可以自己看文档,但是我们只需要学会索引及比对就好了。有点类似于window下面的软件有一个个菜单栏一样,需要自己的鼠标点击来实现一个个功能,在linux下面就是把命令准备好,然后运行。

4.把fastq文件比对到参考基因组

for i in $(seq 1 6) ;do (nohup ~/biosoft/bwa/bwa-0.7.15/bwa  mem -t 5 -M ~/reference/index/bwa/hg19  KPGP-00001_L${i}_R1.fq.gz KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log &);done

这个命令就一句话,但是里面的信息量非常大, 需要熟练掌握linux命令以及shell脚本的语法,但是解析起来也很简单,就是因为我们的fastq文件命名是有规律的,根据规律我构造出一个循环命令,里面的i这个变量会自动扩展成1,2,3,4,5,6依次来用bwa  mem 模式来比对,因为是PE150测序,所以选择这个模式,-M就是选择我们上一步构建好的参考基因组,最后面的 1> 和2>是把软件运行结果输出来,分别是标准输出和标准错误输出,大家可以自行搜索。如果fastq文件的命名发生变化,这个shell脚本是运行不了的,需要临时构建,自己得掌握脚本编写,不然就一个个的比对,手动。

大家可以去看【直播】我的基因组(七):从整体理解全基因组测序数据的变异位点,来了解这个命令的运行结果。

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

1

十二 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

26

【直播】我的基因组(三):抽血送样测序

欢迎收看第三期直播:抽血送样测序。我们从以下几个问题进行本次直播?如果你有什么问题或者不同看法可在微信中留言与作者进行讨论!

FAQ
什么是全基因组测序?

为什么选择二代测序不选择三代测序呢?

为什么要取血液样品测序?

采血过程有什么注意事项?

全基因组测序应该达到何种要求?

因为不了解观看直播的受众知识背景如何,但肯定不会都是生物信息学工程师,所以我还是先把我自己对全基因组测序的了解用通俗的语言描述一下。

    大家都知道,通常正常人都有22对常染色体加上X,Y这两条性染色体。基因组是指生物体所携带的一套完整的单倍体序列,包括全套基因和间隔序列,也就是说我们所描述人类全基因组的时候指的一般是22条常染色体加上X或X、Y的性染色体,它们都是由A、T、C、G碱基组成,总共长度大约是30亿个碱基。全基因组测序,就是检测出全部30亿个碱基对是如何排列的,从第一个到第30亿个,一个都不落下。

我们可以直接在NCBI或者其它网站下载人类的参考基因组序列,因为基因组的探索是持续不断的进化过程,到目前为止,用的是hg38版本的参考基因组。如果你下载了hg38.fa这样的参考基因组记录文件,用文本编辑器打开,就会发现就是ATCG这样的字符。而如果我们想测自己的全基因组,就必须先了解人类的参考全基因组是如何得到的(自行搜索即可)。

到目前为止,主流测序技术仍然是illumina公司的二代测序技术,成本低廉,虽然测序长度就150~250个碱基,对A、T、C、G含量不平衡区域几乎无解,但足以应对大部分的常规分析。我们的全基因组测序一般就用二代测序技术即可,取几百万的细胞破碎后,把所有的染色体随机打断成小片段,一个个的测序,测上亿个片段得到的数据量就很可观了。不过做测序的仪器会比较昂贵,现在最经济,通量最大的是illumina的X10,国内有这个仪器的公司不多!

而三代测序是单分子测序,可以不用进行PCR,就直接对每一条DNA分子的单独测序,但在市场应用中并不是特别成熟,长度可以达到10,000~50,000 nt,但准确度比较低,其次成本也很高。因为读长长,三代测序在组装方面有着很大的优势。

相信大家很多都听说过只需要吐一口唾沫就可以进行基因检测,取一个毛发就可以进行亲子鉴定,但是我们要做的是全基因组测序,所以抽静脉外周血是最好的选择。我这里简单的说说我自己的理解,为什么选择血液而不选择唾沫和毛发做全基因组检测。

唾液,从唾液中提取DNA肯定是没有问题的,不少科研小组都在改进相应的实验操作方法,也很容易搜索到各种论文及新闻报道。然而提取得到的DNA中必然会有口腔微生物的DNA混杂,虽然我们可以把它们一起测序,然后通过比对到人类的参考基因组的手段来去除这些污染,但是我们的测序数据量是有限的,哪怕是口腔微生物的DNA测了序也是需要付费的,大致会有3成的数据被浪费掉。

毛发,从毛发里面提取DNA的技术本身就不成熟,而且提取效率低,再说了也同样有微生物的污染问题。

所以只有血液样品,可以完全避免微生物污染的问题,除非受试者正患有菌血症,否则血液里面是不会有微生物的,而且从血液里面分离白细胞然后提取DNA的技术是非常成熟的,具体实验操作守则见文末。

既然已经确定血液作为样品进行全基因组测序,那么我们就只需要搞定这一个步骤——抽血。每年一次的体检,还有其它时候做过的血常规,我们都有抽过血,很简单,用橡皮筋把胳膊系住,暴露出静脉,然后插入针头,连接采血管,血压自然把学挤到针管尽头的采血管即可(其实我也只是嘴上说说,我自己肯定做不了,也不可能找朋友帮忙,这必须是技术活呀)。

本来想着医院护士做起来那么简单,就直接去医院挂号好了。但连跑了3家医院都被告知不合规,要我们断了这个念头,医院所抽血液属于医疗用品,个人是无权带走的,更别提送去测序公司测序了,当然除非去那种不合规的私人小门诊。当时30斤的干冰都运到了医院,就准备抽完血直接把采血管放到干冰中直接顺丰送到公司测序的,这个壁碰的可真疼!幸运的是,后来在朋友的帮助下找到了主治医师直接安排了一个护士给我用采血管抽血了,整个过程耗时不到5分钟(再次向帮助我的朋友们表达真挚的感激之情)。

如果有也想自己测全基因组的朋友看到本文,请尤其注意采血送样这个步骤,不是想象中的那么简单,一定要先跟测序公司咨询好流程。而且采血管的种类还不是一般的多,这里需要的是EDTA抗凝管紫色头盖,适用于一般的血液学检验,当然也包括提取DNA来进行测序咯。

那么全基因组重测序应该达到什么样的数据要求,搜索一下就可以看到大部分公司宣传的全基因组测序都是30X,就是平均下来能把我们的30亿个碱基每个都测到30次,因为测序是随机的,必然有一些测序深度高一点,有些低一点。至于为什么选择30X这个标准呢,应该是有一篇文章做过梯度模拟,看看5~60X直接,对遗传变异的发现能力的增长情况如何,就是所谓的饱和度分析,而我们全基因组重测序的分析要点,就是挖掘跟参考基因组不一样的地方,而测序深度的增长伴随的就是成本的增长,根据文献(Sequencing depth and coverage: key considerations in genomic analyses - Nature Reviews (2014))及illumina的解释(Sequencing Coverage Calculation Methods for Human Whole-Genome Sequencing;Calling Sequencing SNPs)表明“平均深度达到30X的时候,可以覆盖基因组的95%”、“ This will lead to confident SNP scores and tolerates areas with somewhat lower coverage ”(有兴趣请回复“文献”查看全文)。所以,对我们来说,30X是最佳的选择,可以以最优的成本来挖掘到足够的遗传变异。但是测序仪产出的数据是有质量好坏的,所以还需要跟测序公司约定测序数据里面质量标准,一般用Q20,Q30的百分比这样的指标来表示。而且一般成熟的测序公司都会有相应的数据分析工程师,有成型的数据处理流程,根据支付费用的多少可以选择不同的服务。

附:外周血白细胞基因组的DNA的提取

 

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

1