十二 11

用limma包对芯片数据做差异分析

下载该R语言包,然后看说明书,需要自己做好三个数据(表达矩阵,分组矩阵,差异比较矩阵),总共三个步骤(lmFit,eBayes,topTable)就可以啦

image001

首先做第一个数据,基因表达矩阵!

自己在NCBI里面可以查到下载地址,然后用R语言读取即可

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

image002

然后做好分组矩阵,如下

image003

然后做好,差异比较矩阵,就是说明你想把那些组拿起来做差异分析,如下

image004

最后输出结果:

我进行了6次比较,所以会输出6次比较结果

image005

最后打开差异结果,解读,说明书如下!

 

 

image006

在我的github有完整代码

 

十二 08

下载最新版的KEGG信息,并且解析好

打开官网:http://www.genome.jp/kegg-bin/get_htext?hsa00001+3101

http://www.genome.jp/kegg-bin/get_htext#A1 (这个好像打不开)

可以在里面找到下载链接

image001

下载得到文本文件,可以看到里面的结构层次非常清楚,

image002

C开头的就是kegg的pathway的ID所在行,D开头的就是属于它的kegg的所有的基因

A,B是kegg的分类,总共是6个大类,42个小类

grep ^A hsa00001.keg

A<b>Metabolism</b>

A<b>Genetic Information Processing</b>

A<b>Environmental Information Processing</b>

A<b>Cellular Processes</b>

A<b>Organismal Systems</b>

A<b>Human Diseases</b>

也可以看到,到目前为止(2015年12月8日20:26:57),共有343个kegg的pathway信息啦

image003

接下来我们就把这个信息解析一下:

perl -alne '{if(/^C/){/PATH:hsa(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' hsa00001.keg >kegg2gene.txt

这样就得到了

image004

但是我发现了一个问题,有些通路竟然是没有基因的,我不是很明白为什么?

C    04030 G protein-coupled receptors [BR:hsa04030]

C    01020 Enzyme-linked receptors [BR:hsa01020]

C    04050 Cytokine receptors [BR:hsa04050]

C    03310 Nuclear receptors [BR:hsa03310]

C    04040 Ion channels [BR:hsa04040]

C    04031 GTP-binding proteins [BR:hsa04031]

那我们来看看kegg数据库更新的情况吧。

首先我们看org.Hs.eg.db这个R包里面自带的数据

Date for KEGG data: 2011-Mar15

org.Hs.egPATH has 5869 entrez genes and 229 pathways

2015年八月我用的时候是 6901 entrez genes and 295 pathways

现在是299个通路,6992个基因

所以这个更新其实很缓慢的,所以大家还在用DAVID这种网络工具做kegg的富集分析结果也差不大!

 

 

详细信息见http://www.genome.jp/kegg/pathway.html

更新信息见:http://www.genome.jp/kegg/docs/upd_map.html

十二 01

我的生物信息学视频上线啦

虽然优酷比较坑,但是他的受众多,大家可以去http://i.youku.com/trainee 上面找到我的全部视频,免费观看!

视频列表是http://www.bio-info-trainee.com/tmp/tutorial/video_list.html ,能顺利点击的链接代表视频录制完毕。

这次视频课程的大纲是http://www.bio-info-trainee.com/tmp/tutorial/syllabus.htm,但是很有可能会修改!

声明
大家好,欢迎观看由生信菜鸟团举办的生物信息学公益视频课程,我是主讲人Jimmy。
本课程仅面向那些生物出身却想转生物信息学数据分析的同学,其它人均不在考虑范围。
本课程在中国大陆录制,会尽量遵循中国大陆法律,如有法律问题,请联系我的律师,谢谢!
我会尽量保证视频中知识点的准确无误,如有任何讲解不当之处,欢迎批评指正!
发邮件联系我(jmzeng1314@outlook.com)即可,没什么必要不要QQ或者微信找我。

还有,本人不差钱,所以视频均免费的,录制完毕后我会放出百度云共享,现在还在测试阶段。

请不要无缘无故的批评我,我不会服务任何人,所以不可能有人是我的顾客,我不需要对你好,爱看不看!

十一 15

hapmap计划资源收集

官网是:http://hapmap.ncbi.nlm.nih.gov/index.html.en

所有的数据都放在ncbi上面:ftp://ftp.ncbi.nlm.nih.gov/hapmap/
现在一般用这个计划的数据主要是拿自己得到的突变数据来跟这个hapmap计划的人种突变数据对比。
有芯片数据,也有WES和WGS数据,随着时间的推进,平台也在更新:
Jul 07 2009 00:00    Directory affy100k
Mar 05 2010 00:00    Directory affy500k
Jun 02 2010 00:00    Directory hapmap3_affy6.0
当然,数据也在更新
Jul 07 2009 00:00    Directory 2005-03_phaseI
Dec 03 2009 00:00    Directory 2005-11_phaseII
Jul 07 2009 00:00    Directory 2007-03
Jul 07 2009 00:00    Directory 2008-03
Jul 07 2009 00:00    Directory 2008-07_phaseIII
Jul 07 2009 00:00    Directory 2008-10_phaseII
Jul 07 2009 00:00    Directory 2009-01_phaseIII
Jul 07 2009 00:00    Directory 2009-02_phaseII+III
Aug 18 2010 00:00    Directory 2010-05_phaseIII
Sep 19 2010 00:00    Directory 2010-08_phaseII+III
数据都被整合好了:
  • Bulk data
    • Genotypes: Individual genotype data submitted to the DCC to date. Phase 3 data is available in PLINK format and HapMap format.
    • Frequencies: Allele & genotype frequencies compiled from genotyping data submitted to the DCC to date. These have also been submitted to dbSNP and should be available in the next dbSNP build.
    • LD Data: Linkage disequilibrium properties D', LOD , R2 compiled from the genotype data to date
    • Phasing Data: Phasing data generated using the PHASE software, compiled from the genotype data to date.
    • Allocated SNPs: dbSNP reference SNP clusters that have been picked and prioritized for genotyping according to several criteria (see info on how SNPs were selected). The file 00README contains per-chromosome SNP counts and further details.
    • CNV Genotypes: CNV data from HapMap3 samples.
    • Recombination rates and Hotspots: Recombination rates and hotspots compiled from the genotyping data.
    • SNP assays: Details about assays submitted to the DCC to date. PCR primers, extension probes etc., specific to each genotyping platform.
    • Perlegen amplicons: Details for mapping Perlegen amplicons to HapMap assayLSID. For primer sequences, see Perlegen's Long Range PCR Amplicon data.
    • Raw data: Raw signal intensity data from HapMap genotypes. Currently includes data from Affymetrix GeneChip 100k and 500k Mapping Arrays.
    • Inferred genotypes: Genotypes inferred using the method of Burdick et al. Nat Genet 38:1002-4.
    • Mitochondrial and chrY haplogroups: Classification of phase I HapMap samples into mtDNA and chrY haplogroups. The distribution shown in Table 4 of the HapMap phase I paper (Nat Genet 38:1002-4) corresponds to unrelated parents in each one of the populations analyzed.
同时也发了很多篇文章:
  • The International HapMap Consortium. Integrating common and rare genetic variation in diverse human populations.
    Nature 467, 52-58. 2010. [Abstract] [PDF] [Supplementary information]
  • The International HapMap Consortium. A second generation human haplotype map of over 3.1 million SNPs.
    Nature 449, 851-861. 2007. [Abstract] [PDF] [Supplementary information]
  • The International HapMap Consortium. A Haplotype Map of the Human Genome. 
    Nature 437, 1299-1320. 2005. [Abstract] [PDF] [Supplementary information]
  • The International HapMap Consortium. The International HapMap Project.
    Nature 426, 789-796. 2003. [Abstract] [PDF] [Supplementary information]
  • The International HapMap Consortium. Integrating Ethics and Science in the International HapMap Project. 
    Nature Reviews Genetics 5, 467 -475. 2004. [Abstract] [PDF]
  • Thorisson, G.A., Smith, A.V., Krishnan, L., and Stein, L.D. The International HapMap Project Web site.
    Genome Research,15:1591-1593. 2005. [Abstract] [PDF]
十一 05

点突变详解

DNA分子中某一个碱基为另一种碱基置换,导致DNA碱基序列异常,是基因突变的一种类型。可分为转换和颠换两类。转换(transitions)是同类碱基的置换(AT→GCGC→AT,颠换(transversions) 是不同类碱基的置换(AT→TACG,GC→CGTA

DNA substitution mutations are of two types. Transitions are interchanges of two-ring purines (A  G) or of one-ring pyrimidines (C  T): they therefore involve bases of similar shape. Transversions are interchanges of purine for pyrimidine bases, which therefore involve exchange of one-ring and two-ring structures.

我们在分析driver mutation的时候会区分各种点突变:

  • 1. CpG transitions
  • 2. CpG transversions
  • 3. C:G transitions
  • 4. C:G transversions
  • 5. A:T transitions
  • 6. A:T transversions

那么,我们有64种密码子,每种密码子都会有9种突变可能,我们如何得到一个所有的突变可能的分类并且打分表格呢?

类似于下面这样的表格:共576行!!!

head category.acgt
AAA>AAT 2 A T 6
AAA>AAC 2 A C 6
AAA>AAG 2 A G 5
AAA>ATA 1 A T 6
AAA>ACA 1 A C 6
AAA>AGA 1 A G 5
AAA>TAA 0 A T 6
AAA>CAA 0 A C 6
AAA>GAA 0 A G 5
AAT>AAA 2 T A 6

tail category.acgt
GGC>GGG 2 C G 2
GGG>AGG 0 G A 3
GGG>TGG 0 G T 4
GGG>CGG 0 G C 4
GGG>GAG 1 G A 3
GGG>GTG 1 G T 4
GGG>GCG 1 G C 4
GGG>GGA 2 G A 3
GGG>GGT 2 G T 4
GGG>GGC 2 G C 4

我本来以为这是一件很简单的事情,写起来,才发现好麻烦

1Capture

 

里面用到的一个函数如下:就是判断突变属于上述六种的哪一种!

2

参考:https://www.mun.ca/biology/scarr/Transitions_vs_Transversions.html

https://en.wikipedia.org/wiki/Transversion

http://www.uvm.edu/~cgep/Education/Mutations.html

突变,也称作单碱基替换(single base substitution),指由单个碱基改变发生的突变

可以分为转换(transitions)和颠换(transversions)两类。

转换:嘌呤和嘌呤之间的替换,或嘧啶和嘧啶之间的替换。

颠换:嘌呤和嘧啶之间的替换。

方便理解下面再附上一张示意图,如下:

 

Transitions_vs_Transversions

 

十一 05

使用mutsig软件来找驱动基因

从数以万计的突变里面找到driver mutation这个课题很大,里面的软件我接触的就有十几个了,但是我尝试了其中几个,总是无法运行成功,不知道为什么,终于今天成功了一个,就是mutsig软件! 其实关于突变数据找driver mutation ,台湾一个大学做了一个数据库DriverDB http://ngs.ym.edu.tw/driverdb/: 还因此发了一篇文章:http://nar.oxfordjournals.org/content/early/2013/11/07/nar.gkt1025.full.pdf,挺不错的!

关于driver mutation的理论最近也进化了很多,算是比较完善了吧,但是我一直没时间静下心来好好补充理论知识,很多软件,都只是用过,很多数据,也只是处理了一下,不知道为什么要去做,╮(╯▽╰)╭扯远了,开始谈这个软件吧!

mutsig软件是broadinstitute出品的,所以可靠性非常好咯,来源于一篇nature文章:http://www.nature.com/nature/journal/v505/n7484/full/nature12912.html,而该软件的地址是:http://www.broadinstitute.org/cancer/cga/mutsig_run 需要简单注册才能下载的。

该nature文章是这样描述这个软件的优点的:We used the most recent version of the MutSig suite of tools, which looks for three independent signals: highmutational burden relative to background expectation, accounting for heterogeneity; clustering of mutations within the gene; and enrichment of mutations in evolutionarily conserved sites. Wecombined the significance levels (P values) fromeach test to obtain a single significance level per gene (Methods).

这个软件需要安装matlab环境才能使用,所以我前面就写了教程,如何安装!http://www.bio-info-trainee.com/?p=1166

如果已经安装好了matlab环境,那么直接下载这个软件就可以使用了,软件解压就OK拉,而且人家还提供了测试文件!

Capture4

软件下载后,解压可以看到里面的一个脚本,软件说明书写的非常简单,当然,使用这个软件也的确非常简单:

run_MutSigCV.sh <path_to_MCR> mutations.maf coverage.txt covariates.txt output.txt 即可,其中所有的数据都是可以下载的,

运行完了测试数据, 就证明你的软件安装没有问题啦!如果你只有突变数据的maf格式,maf格式可以参考:https://www.biostars.org/p/69222/ ,也可以使用该软件:如下

run_MutSigCV.sh <path_to_MCR> my_mutations.maf exome_full192.coverage.txt gene.covariates.txt my_results mutation_type_dictionary_file.txt chr_files_hg19

Capture5

上面三个zip文件,都是可以在mutsig软件官网找到下载链接的,是必须下载的!使用很简单,就一个命令即可,但是把你的vcf突变数据做成该软件需要的maf格式,是一个难题!

十一 05

在linux系统里面安装matlab运行环境mcr

matlab毕竟是收费软件,而且是有界面的。所以搞生物信息的都用R和linux替代了,但是很多高大上的单位,比如大名鼎鼎的broadinstitute,是用matlab的,所以他们开发的程序也会以matlab代码的形式发布。但是考虑到大多研究者用不起matlab,或者不会用,所以就用linux系统里面安装matlab运行环境来解决这个问题,我们仍然可以把人家写的matlab程序,在linux命令行下面,当做一个脚本来运行!

比如,这次我就需要用broadinstitute的一个软件:Mutsig,找cancer driver gene的,http://www.broadinstitute.org/cancer/cga/mutsig_run,但是我看了说明才发现,它是用matlab写的,所以我要想在我的服务器用,就必须按照安装matlab运行环境,在官网可以下载:http://www.mathworks.com/products/compiler/mcr/

Capture3

我这里选择的是R2013a (8.1),下载之后解压是这样的,压缩包约四百多M

Capture1

然后直接在解压后的目录里面运行那个install即可,然后如果你的linux可以传送图像,那么就会想安装windows软件一样方便!如果你的linux是纯粹的命令行,那么,就需要一步步的命令行交互,选择安装地址,等等来安装了。

记住你安装之后,会显示一些环境变量给你,请千万要记住,然后自己去修改自己的环境变量,如果你忘记了,就需要搜索来解决环境变量的问题啦!安装之后是这样的:

Capture2

请记住你的安装目录,以后你运行其它matlab相关的程序,都需要把这个安装目录,当做一个参数传给你的其它程序的!!!

如果你没有设置环境变量,就会出各种各样的错误,用下面这个脚本可以设置

其中MCRROOT一般是$path/biosoft/matlab_running/v81/ 这样的东西,请务必注意,LD_LIBRARY_PATH非常重要,非常重要,非常重要!!!!

MCRROOT=$1
echo ---
LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64;
MCRJRE=${MCRROOT}/sys/java/jre/glnxa64/jre/lib/amd64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/native_threads ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/server ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/client ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE} ;
XAPPLRESDIR=${MCRROOT}/X11/app-defaults ;
export LD_LIBRARY_PATH;
export XAPPLRESDIR;
echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH};

如果你没有设置正确,那么会报一下的错误!

error while loading shared libraries: libmwmclmcrrt.so.8.1: cannot o pen shared object file: No such file or directory

error while loading shared libraries: libmwlaunchermain.so: cannot o pen shared object file: No such file or directory

等等!!!

 

十一 01

WES(七)看de novo变异情况

de novo变异寻找这个也属于snp-calling的一部分,但是有点不同的就是该软件考虑了一家三口的测序文件,找de novo突变

软件有介绍这个功能:http://varscan.sourceforge.net/trio-calling-de-novo-mutations.html

而且还专门有一篇文章讲ASD和autism与de novo变异的关系,但是文章不清不楚的,没什么意思

Trio Calling for de novo Mutations

image001

Min coverage:   10

Min reads2:     4

Min var freq:   0.2

Min avg qual:   15

P-value thresh: 0.05

Adj. min reads2:        2

Adj. var freq:  0.05

Adj. p-value:   0.15

Reading input from trio.filter.mpileup

1371416525 bases in pileup file (137M的序列)

83123183 met the coverage requirement of 10 (其中有83M的测序深度大于10X)

145104 variant positions (132268 SNP, 12836 indel) (共发现15.5万的变异位点)

4403 were failed by the strand-filter

139153 variant positions reported (126762 SNP, 12391 indel)

502 de novo mutations reported (376 SNP, 126 indel) (真正属于 de novo mutations只有502个)

1734 initially DeNovo were re-called Germline

12 initially DeNovo were re-called MIE

3 initially DeNovo were re-called MultAlleles

522 initially MIE were re-called Germline

1 initially MIE were re-called MultAlleles

3851 initially Untransmitted were re-called Germline

然后我看了看输出的文件trio.mpileup.output.snp.vcf

软件是这样解释的:The output of the trio subcommand is a single VCF in which all variants are classified as germline (transmitted or untransmitted), de novo, or MIE.

  • FILTER - mendelError if MIE, otherwise PASS
  • STATUS - 1=untransmitted, 2=transmitted, 3=denovo, 4=MIE
  • DENOVO - if present, indicates a high-confidence de novo mutation call

里面的信息量好还是不清楚

我首先对我们拿到的trio.de_novo.mutaion.snp.vcf文件进行简化,只看基因型!
head status.txt   (顺序是dad,mom,child)
STATUS=2 0/0 0/1 0/1
STATUS=2 1/1 1/1 1/1
STATUS=2 0/1 0/0 0/1
STATUS=2 1/1 1/1 1/1
STATUS=1 0/1 0/0 0/0
STATUS=1 0/1 0/0 0/0
STATUS=2 1/1 1/1 1/1
STATUS=2 1/1 1/1 1/1
STATUS=2 1/1 1/1 1/1
STATUS=2 0/1 0/1 0/1
那么总结如下:
  26564 STATUS=1 无所以无 (0/0 0/1 0/0或者 0/1 0/0 0/0等等)
  97764 STATUS=2 有所以有 (1/1 1/1 1/1 或者0/1 0/1 1/1等等)
    385 STATUS=3 无中生有 (0/0 0/0 0/1 或者0/0 0/0 1/1)
   1485 STATUS=4 有中生无 (1/1 0/1 0/0 等等)

我用annovar注释了一下

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  trio.mpileup.output.snp.vcf > trio.snp.annovar

/home/jmzeng/bio-soft/annovar/annotate_variation.pl -buildver hg19  --geneanno --outfile  trio.snp.anno trio.snp.annovar /home/jmzeng/bio-soft/annovar/humandb

结果是:

A total of 132268 locus in VCF file passed QC threshold, representing 132809 SNPs (86633 transitions and 46176 transversions) and 3 indels/substitutions

可以看到最后被注释到外显子上面的突变有两万多个

23794  284345 3123333 trio.snp.anno.exonic_variant_function

这个应该是非常有意义的,但是我还没学会后面的分析。只能先做到这里了

 

十一 01

WES(六)用annovar注释

使用annovar软件参考自:http://www.bio-info-trainee.com/?p=641

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample3.varscan.snp.vcf > Sample3.annovar

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample4.varscan.snp.vcf > Sample4.annovar

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample5.varscan.snp.vcf > Sample5.annovar

然后用下面这个脚本批量注释

image001

Reading gene annotation from /home/jmzeng/bio-soft/annovar/humandb/hg19_refGene.txt ... Done with 50914 transcripts (including 11516 without coding sequence annotation) for 26271 unique genes

最后查看结果可知,真正在外显子上面的突变并不多

23515 Sample3.anno.exonic_variant_function

23913 Sample4.anno.exonic_variant_function

24009 Sample5.anno.exonic_variant_function

annovar软件就是把我们得到的十万多个snp分类了,看看这些snp分别是基因的哪些位置,是否引起蛋白突变

downstream

exonic

exonic;splicing

intergenic

intronic

ncRNA_exonic

ncRNA_intronic

ncRNA_splicing

ncRNA_UTR3

ncRNA_UTR5

splicing

upstream

upstream;downstream

UTR3

UTR5

UTR5;UTR3

 

十一 01

WES(五)不同软件比较

主要是画韦恩图看看,参考:http://www.bio-info-trainee.com/?p=893

对合并而且过滤的高质量snp信息来看看四种不同的snp calling软件的差异

我们用R语言来画韦恩图

image001

可以看出不同软件的差异还是蛮大的,所以我只选四个软件的公共snp来进行分析

首先是sample3

image002

然后是sample4

image003

然后是sample5

image004

可以看出,不同的软件差异还是蛮大的,所以我重新比较了一下,这次只比较,它们不同的软件在exon位点上面的snp的差异,毕竟,我们这次是外显子测序,重点应该是外显子snp

image005

然后我们用同样的程序,画韦恩图,这次能明显看出来了,大部分的snp位点都至少有两到三个软件支持

所以,只有测序深度达到一定级别,用什么软件来做snp-calling其实影响并不大。

 

image006

 

十一 01

WES(四)不同个体的比较

3-4-5分别就是孩子、父亲、母亲

我对每个个体取他们的四种软件的公共snp来进行分析,并且只分析基因型,看看是否符合孟德尔遗传定律

image001

结果如下:

image002

粗略看起来好像很少不符合孟德尔遗传定律耶

然后我写了程序计算

image003

总共127138个可以计算的位点,共有18063个位点不符合孟德尔遗传定律,而且它们在染色体的分布情况如下

我检查了一下,不符合的原因,发现我把

chr1 100617887 C T:DP4=0,0,36,3 T:1/1:40 T:1/1:0,40:40 miss T:DP4=0,0,49,9 T:1/1:59 T:1/1:0,58:59 miss T:DP4=0,0,43,8 T:1/1:53 T:1/1:0,53:53 T:1/1:50

计算成了chr1 100617887 C 0/0 0/0 1/1 所以认为不符合,因为我认为只有四个软件都认为是snp的我才当作是snp的基因型,否则都是0/0

那么我就改写了程序,全部用gatk结果来计算。这次可以计算的snp有个176036,不符合的有20309,而且我看了不符合的snp的染色体分布,Y染色体有点异常

image004 image005

但是很失败,没什么发现!

 

 

十一 01

WES(三)snp-filter

其中freebayes,bcftools,gatk都是把所有的snp细节都call出来了,可以看到下面这些软件的结果有的高达一百多万个snp,而一般文献都说外显子组测序可鉴定约8万个变异!

image001

这样得到突变太多了,所以需要过滤。这里过滤的统一标准都是qual大于20,测序深度大于10。过滤之后的snp数量如下

image002

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample3.bcftools.vcf >Sample3.bcftools.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample4.bcftools.vcf >Sample4.bcftools.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample5.bcftools.vcf >Sample5.bcftools.vcf.filter

 

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample3.freebayes.vcf > Sample3.freebayes.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample4.freebayes.vcf > Sample4.freebayes.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample5.freebayes.vcf > Sample5.freebayes.vcf.filter

 

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample3.gatk.UG.vcf  >Sample3.gatk.UG.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample4.gatk.UG.vcf  >Sample4.gatk.UG.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample5.gatk.UG.vcf  >Sample5.gatk.UG.vcf.filter

 

perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample3.varscan.snp.vcf >Sample3.varscan.snp.vcf.filter

perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample4.varscan.snp.vcf >Sample4.varscan.snp.vcf.filter

perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample5.varscan.snp.vcf >Sample5.varscan.snp.vcf.filter

这样不同工具产生的snp记录数就比较整齐了,我们先比较四种不同工具的call snp的情况,然后再比较三个人的区别。

然后写了一个程序把所有的snp合并起来比较

image003

得到了一个很有趣的表格,我放在excel里面看了看 ,主要是要看生物学意义,但是我的生物学知识好多都忘了,得重新学习了 image004

十一 01

WES(二)snp-calling

准备文件:下载必备的软件和参考基因组数据

1、软件

ps:还有samtools,freebayes和varscan软件,我以前下载过,这次就没有再弄了,但是下面会用到

2、参考基因组

3、参考 突变数据

第一步,下载数据

第二步,bwa比对

第三步,sam转为bam,并sort好

第四步,标记PCR重复,并去除

第五步,产生需要重排的坐标记录

第六步,根据重排记录文件把比对结果重新比对

第七步,把最终的bam文件转为mpileup文件

第八步,用bcftools 来call snp

第九步,用freebayes来call snp

第十步,用gatk     来call snp

第十一步,用varscan来call snp

下面的图片是按照顺序来的,我就不整理了

image001 image002 image003 image004 image005

image006 image007 image008 image009 image010 image011 image012 image013 image014 image015 image016 image017 image018 image019 image020 image021

 

十一 01

基因及外显子到底占基因组多少比例

很多文献都提到外显子组的序列仅占全基因组序列的1%左右,但大多数与疾病相关的变异位于外显子区。通过外显子组测序可鉴定约8万个变异,全基因组测序可鉴定300万个变异,因此与全基因组测序相比,外显子组测序不仅费用较低,数据阐释也更为简单。外显子组测序技术以其经济有效的优势广泛应用于孟德尔遗传病、罕见综合征及复杂疾病的研究,并于2010年被Science杂志评为十大突破之一。所以我一直想验证一下外显子到底占基因组什么样的百分比,是哪些位点。

首先我们拿到外显子记录和基因信息,下面是通过R来从bioconductor中心提取数据

library("TxDb.Hsapiens.UCSC.hg19.knownGene")

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

txdb

exon_txdb=exons(txdb)

genes_txdb=genes(txdb)

tmp=as.data.frame(exon_txdb)

write.table(tmp,'exon.pos',row.names=F)

tmp=as.data.frame(genes_txdb)

write.table(tmp,'gene.pos',row.names=F)

首先我们看看我刚才从TxDb.Hsapiens.UCSC.hg19.knownGene包里面提取的外显子记录文件exon.pos

image001

我简单用脚本统计了一下:perl -alne '{$tmp+=$F[3]} END{print $tmp}' exon.pos

共289970个外显子记录,长度共98759283bp,也就是98.5Mb的区域都是外显子,感觉有点不大对,毕竟真正的外显子也就三十多M的,所以必须有些外显子是并列关系,也就是有的外显子包含着另外一些外显子,然后我写脚本检查了一下,的确太多了,这样的重复!

image002

也可以去NCBI里面拿到 consensus coding sequence (CCDS)记录:ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/

那么我再看看NCBI里面拿到 consensus coding sequence (CCDS)记录CCDS.20150512.txt文件。

共33140行,一个基因一行,所以需要格式化成外显子文件,简单看了一下文件的列,很容易格式化。每一行的信息如下

1 NC_000001.11 SAMD11 148398 CCDS2.2 Public + 925941 944152 [925941-926012, 930154-930335, 931038-931088, 935771-935895, 939039-939128, 939274-939459, 941143-941305, 942135-942250, 942409-942487, 942558-943057, 943252-943376, 943697-943807, 943907-944152] Identical

用下面这个脚本格式化

image003

格式化后的文件如下

image004

外显子记录增加到了346493个,我再用单行命令计算长度,是56321786bp,这次只有56M了,还是有点大,所以应该是还有重复!

所以,我写了一个脚本来精确计算外显子的总长度,不能那样马马虎虎的用单行命令了。

image005

这次才是34729283bp,也就是约35M,根很多文献里面说的都一样!

同理,我统计了一下外显子的侧翼长度

54692160 外显子加上前后50bp

73066288  外显子加上前后100bp

90362533  外显子加上前后150bp

 

而hg19版本基因组里面有着entrez gene ID号的基因是23056个基因,所以我接下来探究一下这些基因的信息!

我们首先看看基因与基因之间的交叉情况

其中有12454266bp的位点,是多个外显子共有的,可能是一个基因的多个外显子,或者是不同基因的外显子

我们主要看看不同基因的交叉情况

不同基因交叉情况共516种,共涉及498种基因!

ZNF341   NFS1

ZNF397   ZSCAN30

ZNF428   SRRM5

ZNF436-AS1    RPL11

ZNF578   ZNF137P

ZNF619   ZNF621

ZNF816   ZNF816-ZNF321P

ZSCAN26  ZSCAN31

ZSCAN5A  ZNF542P

ZZEF1    ANKFY1

我随便在数据库里面搜索一下验证,发现这些基因的确是交叉重叠的

 

30

R语言也可以用hash啦

本来想着是如何习惯R里面的list对象,来实现hash的功能,无意中发现了居然还有这个包

https://cran.r-project.org/web/packages/hash/hash.pdf

我简单看了看文档,的确跟perl里面的hash功能类似,非常好用。

创建一个hash很简单

h <- hash( keys=letters, values=1:26 )

h <- hash( letters, 1:26 )

类似于下面的列表

L=setNames(as.list(LETTERS),1:26)

如果是列表要提取keys需要用names函数,如果需要提取values,需要用unlist函数,而用了hash之后,这两个函数就可以独自运行啦

还有很多其它函数,其实在list中也可以实现,主要是看你对哪种语法更熟悉,感觉自己现在编程能力差不多算是小有所成了,看什么都一样了,要是以前学R的时候看到了这个hash包,我肯定会很兴奋的!

clear signature(x = "hash"): Remove all key-value pairs from hash
del signature(x = "ANY", hash = "hash"): Remove specified key-value pairs from hash
has.key signature(key = "ANY", hash = "hash"): Test for existence of key
is.empty signature(x = "hash"): Test if no key-values are assigned
length signature(x = "hash"): Return number of key-value pairs from the hash
keys signature(hash = "hash"): Retrieve keys from hash
values signature(x = "hash"): Retrieve values from hash
copy signature(x = "hash"): Make a copy of a hash using a new environment.
format signature(x = "hash"): Internal function for displaying hash

 

30

R里面list对象和AnnDbBimap对象区别

AnnDbBimap对象是R里面的bioconductor系列包的基础对象,在探针数据里面会包装成ProbeAnnDbMap,跟go通路相关又是GOTermsAnnDbBimap对象。

但是它都是AnnDbBimap对象衍生过来的

主要存在于芯片系列的包和org系列的包,其实AnnDbBimap对象就是list对象的包装,比如下面这些例子:

ls("package:hgu133plus2.db")

 [1] "hgu133plus2"              "hgu133plus2.db"           "hgu133plus2_dbconn"
 [4] "hgu133plus2_dbfile"       "hgu133plus2_dbInfo"       "hgu133plus2_dbschema"
 [7] "hgu133plus2ACCNUM"        "hgu133plus2ALIAS2PROBE"   "hgu133plus2CHR"
[10] "hgu133plus2CHRLENGTHS"    "hgu133plus2CHRLOC"        "hgu133plus2CHRLOCEND"
[13] "hgu133plus2ENSEMBL"       "hgu133plus2ENSEMBL2PROBE" "hgu133plus2ENTREZID"
[16] "hgu133plus2ENZYME"        "hgu133plus2ENZYME2PROBE"  "hgu133plus2GENENAME"
[19] "hgu133plus2GO"            "hgu133plus2GO2ALLPROBES"  "hgu133plus2GO2PROBE"
[22] "hgu133plus2MAP"           "hgu133plus2MAPCOUNTS"     "hgu133plus2OMIM"
[25] "hgu133plus2ORGANISM"      "hgu133plus2ORGPKG"        "hgu133plus2PATH"
[28] "hgu133plus2PATH2PROBE"    "hgu133plus2PFAM"          "hgu133plus2PMID"
[31] "hgu133plus2PMID2PROBE"    "hgu133plus2PROSITE"       "hgu133plus2REFSEQ"
[34] "hgu133plus2SYMBOL"        "hgu133plus2UNIGENE"       "hgu133plus2UNIPROT"

那么我们可以随便挑选包中的一个数据分析一下

x <- hgu133plus2SYMBOL

xlist=as.list(x)

我们查看X对象,可知,它是object of class "ProbeAnnDbMap",而这个对象,就是常见的list对象,被包装了一下, 只有我们明白了它和list对象到底有什么区别,才算是真正搞懂了这一系列包

但是这个ProbeAnnDbMap对象,在GO等包里面会被包装成更复杂的对象-GOTermsAnnDbBimap,但是对他们的理解都大同小异。

我们先回顾一下在R语言里面的list的基础知识:

R 的 列表(list)是一个以对象的有序集合构成的对象。 列表中包含的对象又称为它的分量(components)。

分量可以是不同的模式或者类型, 如一个列表可以包括数值向量,逻辑向量,矩阵,复向量, 字符数组,函数等等

如果你会perl的话,可以把它理解成hash。

分量常常会被编号的(numbered),并且可以利用它来 访问分量

列表的分量可以被命名,这种情况下 可以通过名字访问。

特别要注意一下 Lst[[1]] 和 Lst[1] 的差别。 [[...]] 是用来选择单个 元素的操作符,而 [...] 是一个更为一般 的下标操作符。

因此前者得到的是列表 Lst 中的第一个对象, 并且含有分量名字的命名列表(named list)中分量的名字会被排除在外的。

后者得到的则是列表 Lst 中仅仅由第一个元素 构成的子列表。如果是命名列表, 分量名字会传给子列表的。

那么接下来我们就看看x和xlist的区别。

它们里面的数据都是affymetrix公司出品的人类的hgu133plus2芯片的探针ID与基因symbol的对应关系

如果我想拿到所有的探针ID,那么对于AnnDbBimap对象,需要用mappedkeys(x),对于普通的list对象,需要用names(xlist).

对于普通的list对象,如果我想看前几个元素,直接head就可以了,但是对于AnnDbBimap对象,数据被封装好了,需要先as.list,然后才能head

mapped_probes <- mappedkeys(x)

PID2=head(mapped_probes)

那么,如果我们想根据以下探针ID来查看它们在这些数据里面被对应着哪些基因symbol 呢?

> PID2 #这是一串探针ID,后面的操作都需要用的

[1] "1053_at"   "117_at"    "121_at"    "1255_g_at" "1316_at"   "1320_at"

如果是对于AnnDbBimap对象,我们可以用mget函数来操作,取多个探针的对应基因symbol

accnum <- mget(PID2, env=hgu133plus2ACCNUM);

gname <- mget(PID2, env=hgu133plus2GENENAME)

gsymbol <- mget(PID2, env=hgu133plus2SYMBOL)

mget函数返回的就是普通的list函数了,可以直接查看了。

如果是对于普通的list对象,我们想取多个探针的对应基因symbol也是非常简单的,xlist[PID2]即可。

那么我们不禁有问了,既然它们两个功能完全一样,何苦把它包装成一个对象了,我直接操作list对象不就好了,学那么多规矩干嘛?

所以,重点就来了

>  length(mapped_probes)

[1] 42125

> length(names(xlist))

[1] 54675

看懂了吗?

 

但是,事实上用处也不大,你觉得下面这两个有区别吗?

SYMBOL <- AnnotationDbi::as.list(hgu133plus2SYMBOL)

SYMBOL <- SYMBOL[!is.na(SYMBOL)];

 

x <- hgu133plus2SYMBOL

mapped_probes <- mappedkeys(x)

SYMBOL <- AnnotationDbi::as.list(x[mapped_probes])

 

PS,在R里面创建一个list对象是非常简单的

The setNames() built-in function makes it easy to create a hash from given key and value lists.(Thanks to Nick K for the better suggestion.)

Usage: hh <- setNames(as.list(values), keys)

Example:

players <- c("bob", "tom", "tim", "tony", "tiny", "hubert", "herbert")
rankings <- c(0.2027, 0.2187, 0.0378, 0.3334, 0.0161, 0.0555, 0.1357)
league <- setNames(as.list(rankings), players)

29

gene的各种ID转换终结者-bioconductor系列包

经常会有人问这样的问题I have list of 10,000 Entrez IDs and i want to convert the multiple Entrez IDs into the respective gene names. Could someone suggest me the way to do this?等等类似的基因转换,能做的基因转换的方法非常多,以前不懂编程的时候,都是用各种网站,而最常用的就是ensembl的biomart了,它支持的ID非常多,高达几百种,好多ID我到现在都不知道是什么意思。

现在学会编程了,我比较喜欢的是R的一些包,是bioconductor系列,一般来说,其中有biomart,org.Hs.eg.db,annotate,等等。关于biomart我就不再讲了,我前面的博客至少有七八篇都提到了它。本次我们讲讲简单的, 我就以把gene entrez ID转换为gene symbol 为例子把。

当然,首先要安装这些包,并且加载。

if("org.Hs.eg.db" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("org.Hs.eg.db")}
suppressMessages(library(org.Hs.eg.db))  我比较喜欢这样加载包

library(annotate) #一般都是这样加载包

如果是用org.Hs.eg.db包,首先你只需要读取你的待转换ID文件,构造成一个向量,tmp,然后只需要symbols <- org.Hs.egSYMBOL[as.character(tmp)]就可以得到结果了,返回的symbols是一个对象,需要用toTable这个函数变成数据框。但是这样转换容易出一些问题,比如如果你的输入数据tmp,里面含有一些无法转换的gene entrez ID,就会报错。

而且它支持的ID转换很有限,具体看看它的说明书即可:https://www.bioconductor.org/packages/release/data/annotation/manuals/org.Hs.eg.db/man/org.Hs.eg.db.pdf

org.Hs.eg.db
org.Hs.eg_dbconn
org.Hs.egACCNUM
org.Hs.egALIAS2EG
org.Hs.egCHR
org.Hs.egCHRLENGTHS
org.Hs.egCHRLOC
org.Hs.egENSEMBL
org.Hs.egENSEMBLPROT
org.Hs.egENSEMBLTRANS
org.Hs.egENZYME
org.Hs.egGENENAME
org.Hs.egGO
org.Hs.egMAP
org.Hs.egMAPCOUNTS
org.Hs.egOMIM
org.Hs.egORGANISM
org.Hs.egPATH
org.Hs.egPMID
org.Hs.egREFSEQ
org.Hs.egSYMBOL
org.Hs.egUCSCKG
org.Hs.egUNIGENE
org.Hs.egUNIPROT

如果是用annotate包,首先你还是需要读取你的待转换ID文件,构造成一个向量,tmp,然后用getSYMBOL(as.character(tmp), data='org.Hs.eg')这样直接就返回的还是以向量,只是在原来向量的基础上面加上了names属性。说明书:http://www.bioconductor.org/packages/3.3/bioc/manuals/annotate/man/annotate.pdf

然后你可以把转换好的向量写出去,如下:

1 A1BG
2 A2M
3 A2MP1
9 NAT1
10 NAT2
12 SERPINA3
13 AADAC
14 AAMP
15 AANAT
16 AARS

PS:如果是芯片数据,需要把探针的ID转换成gene,那么一般还需要加载特定芯片的数据包才行:

platformDB <- paste(eset.mas5@annotation, ".db", sep="") #这里需要确定你用的是什么芯片
cat("the annotation is ",platformDB,"\n")
if(platformDB %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");tmp=try(biocLite(platformDB))}
library(platformDB, character.only=TRUE)
probeset <- featureNames(eset.mas5)
rowMeans <- rowMeans(exprSet)

library(annotate) # lookUp函数是属于annotate这个包的
EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

29

使用GEOmetadb包来获取对应GEO数据的实验信息

理论上我前面提到的GEOquery包就可以根据一个GSE索引号来获取NCBI提供的所有关于这个GSE索引号的数据了,包括metadata,表达矩阵,soft文件,还有raw data

但是很多时候,那个metadata并不是很整齐,而且一个个下载太麻烦了,所以就需要用R的bioconductor的另一个神奇的包了GEOmetadb

它的示例:http://bioconductor.org/packages/devel/bioc/vignettes/GEOmetadb/inst/doc/GEOmetadb.R

里面还是很多数据库基础知识的
代码托管在github,它的示例代码是这样连接数据库的:
library(GEOmetadb)
if(!file.exists('GEOmetadb.sqlite')) getSQLiteFile()
file.info('GEOmetadb.sqlite')
con <- dbConnect(SQLite(),'GEOmetadb.sqlite')
dbDisconnect(con)
但是一般不会成功,因为这个包把它的GEOmetadb.sqlite文件放在了国外网盘共享,在国内很难访问,推荐大家想办法下载到本地
tmp2
用这个代码就会成功了,需要自己下载GEOmetadb.sqlite文件然后放在指定目录:/path/GEOmetadb.sqlite 需要自己修改
我们的diabetes.GEO.list文件内容如下:
GSE1009
GSE10785
GSE1133
GSE11975
GSE121
GSE12409
那么会产生的表格文件如下:共有32列数据信息,算是蛮全面的了
tmp
28

模拟Y染色体测序判断,并比对到X染色体上面,看同源性

首先下载两条染色体序列

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrX.fa.gz;

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrY.fa.gz;

152M Mar 21  2009 chrX.fa

58M Mar 21  2009 chrY.fa

然后把X染色体构建bwa的索引

bwa index chrX.fa

[bwa_index] Pack FASTA... 1.97 sec

[bwa_index] Construct BWT for the packed sequence...

[BWTIncCreate] textLength=310541120, availableWord=33850812

[BWTIncConstructFromPacked] 10 iterations done. 55838672 characters processed.

[BWTIncConstructFromPacked] 20 iterations done. 103157920 characters processed.

[BWTIncConstructFromPacked] 30 iterations done. 145211344 characters processed.

[BWTIncConstructFromPacked] 40 iterations done. 182584528 characters processed.

[BWTIncConstructFromPacked] 50 iterations done. 215797872 characters processed.

[BWTIncConstructFromPacked] 60 iterations done. 245313968 characters processed.

[BWTIncConstructFromPacked] 70 iterations done. 271543920 characters processed.

[BWTIncConstructFromPacked] 80 iterations done. 294853104 characters processed.

[bwt_gen] Finished constructing BWT in 88 iterations.

[bwa_index] 98.58 seconds elapse.

[bwa_index] Update BWT... 0.96 sec

[bwa_index] Pack forward-only FASTA... 0.91 sec

[bwa_index] Construct SA from BWT and Occ... 33.18 sec

[main] Version: 0.7.8-r455

[main] CMD: /lrlhps/apps/bioinfo/bwa/bwa-0.7.8/bwa index chrX.fa

[main] Real time: 141.623 sec; CPU: 135.605 sec

由于X染色体也就152M,所以很快,两分钟解决战斗!

然后模拟Y染色体的测序判断(PE100insert400

209M Oct 28 10:19 read1.fa

209M Oct 28 10:19 read2.fa

模拟的程序很简单

tmp

 

while(<>){
chomp;
$chrY.=uc $_;
}
$j=0;
open FH_L,">read1.fa";
open FH_R,">read2.fa";
foreach (1..4){
for ($i=600;$i<(length($chrY)-600);$i = $i+50+int(rand(10))){
$up = substr($chrY,$i,100);
$down=substr($chrY,$i+400,100);
next unless $up=~/[ATCG]/;
next unless $down=~/[ATCG]/;
$down=reverse $down;
$down=~tr/ATCG/TAGC/;
$j++;
print FH_L ">read_$j/1\n";
print FH_L "$up\n";
print FH_R ">read_$j/2\n";
print FH_R "$down\n";
}
}
close FH_L;
close FH_R;

然后用bwa mem 来比对

bwa mem -t 12 -M chrX.fa read*.fa >read.sam

用了12个线层,所以也非常快

[main] Version: 0.7.8-r455

[main] CMD: /apps/bioinfo/bwa/bwa-0.7.8/bwa mem -t 12 -M chrX.fa read1.fa read2.fa

[main] Real time: 136.641 sec; CPU: 1525.360 sec

643M Oct 28 10:24 read.sam

然后统计比对结果

samtools view -bS read.sam >read.bam

158M Oct 28 10:26 read.bam

samtools flagstat read.bam

3801483 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

2153410 + 0 mapped (56.65%:-nan%)

3801483 + 0 paired in sequencing

1900666 + 0 read1

1900817 + 0 read2

645876 + 0 properly paired (16.99%:-nan%)

1780930 + 0 with itself and mate mapped

372480 + 0 singletons (9.80%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

我自己看sam文件也发现真的同源性好高呀,总共就模拟了380万reads,就有120万是百分百比对上了。

所以对女性个体来说,测序判断比对到Y染色体是再正常不过的了。如果要判断性别,必须要找那些X,Y差异性区段