30

用谷歌搜索来使用ggplot2做可视化(下)

用谷歌搜索来使用ggplot2做可视化(下)

2017-01-30 jimmy 生信菜鸟团

我知道会有续集,但也没想到续集来得这么快!今天收到了一个生信技能树公众账号铁杆粉丝(我们之间有过9次邮件交流)的求助信,下面我们首先一起帮他解决一下碰到的问题。随后和大家分享一下可以提高搜索效率和准确率的Google搜索技巧。

 

Continue reading

30

如何通过Google来使用ggplot2可视化(上)

如何通过Google来使用ggplot2可视化

2017-01-29 jimmy 生信菜鸟团

今天是大年初二,这篇文章我只想传达一点:

没有什么菜鸟级别的生物信息学数据处理是不能通过Google得到解决方案的,如果有,请换个关键词继续Google!


第一部分

首先用两分钟的时间简单介绍一下R语言

因为这个语言是肉丝儿Ross Ihaka)和萝卜特Robert Gentleman)两个人1992年在S语言的基础上发明出来的开源语言,所以叫做R语言。这两个人是统计学教授出身,所以R语言在统计学方面有着纯正的血统!如果你平时的工作和统计相关,你好意思不会点R语言么?

 

Continue reading

23

给初学者的忠告,不要拿一套垃圾数据入门!

垃圾数据对初学者的伤害真的很可怕!

最近在带一些朋友入门,想起了当年自己入门的各种凄惨惨戚戚!

碱基质量值很差,GC不平衡,还有接头,PCR重复也很多,kmer值也很诡异,时间都耗在QC上面了,结果几个月下来,你一个流程都没搞明白,各种查资料,还是在原地打转。 Continue reading

14

转一个Python的安利文章咯!

来自于我们生信技能树论坛的超级版主bioinfo.dong的好文一篇,比较符合我博客的思想,就友情转发一下:

原文链接见:http://www.biotrainee.com/thread-379-1-1.html 

刚接触生信的同学大都有个困惑,知道生物信息可能需要编程,可是选择什么语言呢?有人会说perl啊,Python啊,R啊,java啊,等等等等。目的不一样,选择也不一样,你可以说语言都没有区别,达到目的就行,当然没问题。可是我们也要知道每种语言都有其独特优势,你可以用perl倒腾出矩阵运算,也可以画出想要的图,可是没有R专业;你也可以用R的正则表达式处理文本,可是perl或者Python做正则会更方便一些。这不是比较帖,只是从一个Python体验者的角度来说一下为什么选择Python。我目前的编程组合是Python+R+Shell Scripting。

这篇文章比较适合编程初学者,常年用perl的老司机们可以随便看一下,虽说perl和Python很像,有了一门的基础,学另一门就容易多了,可是真让一个用了几年perl的人彻底换Python还是比较困难的,主要还是习惯问题。最初做生信的人大都以perl作为常用脚本语言。我也是从perl开始的,当年为了申请出国读Bioinformatics,认真把小骆驼书看了一遍。来美国之后的第一个导师刚好是教perl的,我又跟着学了一次,看完导师推荐的《Unix and Perl to the Rescue》,算是巩固加第二次入门。之后一年基本都是用perl来处理数据。一个偶然的机会,同学说一起学学Python吧,听说很好用,于是就在网上找了个教程把题目刷了一遍。虽说入了门,可是每次项目赶时间的时候第一个想到的还是用perl来解决,所以入门很久也没啥长进,我亲爱的同学因为perl用的太好,虽然知道Python很好用,可始终没法狠心转过来,而我因为本身perl学得也只是半斤八两,纠结了一段时间也就彻底放弃perl了。

先说用了很长时间perl再用Python觉得不习惯的点。

(1)首先是动物园的书,《learning perl》真是入门的典范。再看《Learning Python》,几千页,那么厚,我到现在也没法认真看下去。
(2)另外perl语句比较简洁,几个符号就可以讲清楚的,Python可能需要几行,比如按行读取,perl只要while(<>)就可以,而最初学Python的时候,光这个问题就困扰了很久。再比如perl正则匹配的$1, Python是match.group(1)。perl的简洁伴随的缺点是可读性较差,自己的代码写完了都不想再看,更不要说别人写的。
(3)perl的正则表达式是真的非常厉害,我已经不记得是怎么厉害的了,就只记得Python的re module刚开始接触不太好用,不过现在已经感觉不出区别了。
(4)通常一个Python脚本需要很多modules,不熟悉之前会觉得很痛苦,perl就比较少用到,我总共也没用几次,一方面说明我的perl确实学得不好,另一方面可能也真是不太好用,看到就觉得麻烦。但Python的modules一旦熟悉了会大大提高工作效率。

重点说一下Python的优点。Python作为编程语言真正的优势比如面向对象编程(OOP),可移植/扩展/嵌入,强大的爬虫功能,APP开发,web开发等都不在讨论范围之内,只从最实用的角度做一下说明:
(1)简单,适合作为入门语言。很多时候觉得读Python的代码像是在读简单的英文,或者觉得pseudocode稍微一改就可以在Python里run了。Python还规范了很好的写作格式,该缩进的必须缩进,这样更增强了可读性。同时提高了代码重复利用的可能(很多时候perl代码写完就不想读了,三个月不用再回来已经看不太懂了,Python的就可以留着慢慢用。。。)
(2)Python社区活跃。有问题可以很容易搜索到解决方案。我perl的老师现在也转教Python了,问他为什么,他说perl的community不活跃,用Python是一种趋势
(3)作为开源语言,Python有很多非常好用的包,可以最大程度让我们避免把时间浪费在重复造轮子上。刚接触Python的时候我就觉得这简直是perl和R的整合,之前提过Python的scipy,numpy,pandas,matlibplot等等packages使其同样拥有了很强大的统计画图功能,我曾一度弃用R,用Python做所有的数据处理,数据分析和画图。不过现在又将这些工作交回了R,实验室本身是做统计的,用R显得入流一点:-)
(4)Python的jupyter notebook!!!这个是要强力推荐的!!!以前叫ipython notebook。用过R的都知道R Studio。jupyter notebook就是Python的Studio。以前写perl或者Python是不是这样的流程:写好了,存成.pl或.py格式,在shell里python xxx.py或者perl xxx.pl。运行完发现不好,有bug,打开文件找找bug在哪,再运行,还不行,唉,反反复复,好累。有了jupyter notebook你就可以边写边跑边改程序。有任何不确定的地方,都可以在notebook里直接测试,有任何bug都可以在notebook里直接改。简直方便到爆。现在用Anaconda安装jupyter还附赠很多包,方便又实惠。
(5)学好Python可以转行!!!跳出生物坑,奔向美好的互联网坑。前面提到的爬虫,APP开发,web编程都是很实用的技能。许多互联网公司也会专门招Python程序员,比如Google,比如Youtube,比如Dropbox。。。

我本专业是Bioinformatics,需要上一些计算机和统计的研究生课程,还记得算法课上老师第一节课就问,java和C++都会吧,如果不会的话Python总会吧,都不会的话这门课的作业写不了。就因为觉得自己还算会一点Python,把一次学习java的好机会浪费掉了

暂时就想到这么多。说的未必对。都是自己的体会吧。希望对初学者有用~

05

生物学基础知识~CARM和SWI/SNF复合物

因为最近在研究CHIP-seq测序数据处理,发现有些文章重点并不是数据处理本身,而是对生物学基础知识的掌控以及实验设计,这篇文章我重点推荐一下,我略微翻译了一些,笔记如下:

SWI/SNF(BAF) chromatin remodeling complex  染色质重构复合物,以及被广泛发现在各种癌症患者体内均有突变,这个复合物利用ATP水解释放的能量来驱动核小体运动以及调控染色质的结构。
发现历史: 最新是在yeast里面发现它的突变会影响 mating type SWItching 和 导致sucrose nonfermenting 的表型,所以才简称为SWI/SNF,在哺乳动物体内也被广泛研究现在,它被发现参与了细胞分化、增殖、还有各种DNA修复功能的实现。
组成结构:本质是一个蛋白质复合物,约15个亚基,每个亚基都由一个独立的基因转录翻译而来,每个基因都有专门的文章研究过。
SWI-SNF-geneFamily-SMARCC1-BAF155
生物功能:早在1998年就被发现与癌症相关,在GO数据库里面还可以查到相关资料,所以很容易在各种通路分析结果里面看到它的身影
GO-0016514-SWI-SNF-complex
随着癌症基因组测序的进展,至少有8种染色质重构复合物的亚基有recurrent mutation情况,然后有一句话说得特别好:
This has resulted in interest in indentifying mechanisms by which activity of the SWI/SNF complex is regulated, with the hope that such mechanistic understanding may reveal novel opportunities for therapeutic intervention .
CARM1基因是PRMT系列的一种:protein arginine Methyltransferases一直与基因的转录与翻译相关,主要功能就是使得蛋白质的arginine甲基化,这个基因家族目前有9个基因,只有CARM1命名比较奇葩。
PRMT-gene-family-CARM1
其中RM都是 Arginine Methyltransferase的简称,很容易理解,
CARM1这个酶的底物不仅仅包括参与染色质重构,还有基因转录调控,剪切因子,组蛋白乙酰化还有RNA结合蛋白
chromatin remodeling (  SWI/SNF(BAF) chromatin remodeling complex  )
gene transcription   (histone H3(at R17))
histone acetyltransferases  ( p300/CBP )
splicing factors  ( CA150,SAP49,SmB,U1C  )
RNA-binding proteins   ( PABP1, HuR, HuD)
也有实验证明非常多的癌症种类病人都发现了CARM1表达异常升高,而且刺激乳腺癌的恶化,而且是很多癌症相关恶化蛋白的共刺激因子,比如p53,E2F1,NF-kB,b-catenin,steroid hormone receptors. 
但是CARM1到底是如何在乳腺癌体内发生作用的,其中机制并不完全清楚。
25

自学miRNA-seq分析第三讲~公共测序数据下载

前面已经讲到了该文章的数据已经上传到NCBI的SRA数据中心,所以直接根据索引号下载,然后用SRAtoolkit转出我们想要的fastq测序数据即可。下载的数据一般要进行质量控制,可视化展现一下质量如何,然后根据大题测序质量进行简单过滤。所以需要提前安装一些软件来完成这些任务,包括: sratoolkit /fastx_toolkit /fastqc/bowtie2/hg19/miRBase/SHRiMP

下面是我用新服务器下载安装软件的一些代码记录,因为fastx_toolkit /fastqc我已经安装过,就不列代码了,还有miRBase的下载,我在前面第二讲里面提到过,传送门:自学miRNA-seq分析第二讲~学习资料的搜集 Continue reading

05

用snpEFF对vcf格式的突变数据进行注释

这个软件比较重要,尤其是对做遗传变异相关研究的,很多人做完了snp-calling后喜欢用ANNOVAR来进行注释,但是那个注释还是相对比较简单,只能得到该突变位点在基因的哪个区域,那个基因这样的信息,如果想了解更具体一点,就需要更加功能化的软件了,snpEFF就是其中的佼佼者,而且是java平台软件,非常容易使用!而且它的手册写的非常详细:http://snpeff.sourceforge.net/SnpEff_manual.html

官网是:http://snpeff.sourceforge.net/

1       889455  .       G       A       .       .        ## 假设我们的vcf文件里面记录的突变是这个,那么我们可以用snpEFF进行注释,注释得到的信息非常完全!

信息用|符号分割,所有很容易用脚本提取需要的信息

ANN=A|stop_gained|HIGH|NOC2L|ENSG00000188976|transcript|ENST00000327044|protein_coding|7/19|c.706C>T|p.Gln236*|756/2790|706/2250|236/749||,A|downstream_gene_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000487214|processed_transcript||n.*865C>T|||||351|,A|downstream_gene_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000469563|retained_intron||n.*878C>T|||||4171|,A|non_coding_exon_variant|MODIFIER|NOC2L|ENSG00000188976|transcript|ENST00000477976|retained_intron|5/17|n.2153C>T||||||;LOF=(NOC2L|ENSG00000188976|6|0.17);NMD=(NOC2L|ENSG00000188976|6|0.17)

包括突变类型是:non_coding_exon_variant

突变在各种转录本上面,在每个转录本的第几个碱基呀,哪个氨基酸的改变呀,氨基酸第几位呀!

标准突变表示形式是:

突变发生在NOC2L这个基因上面,它的ensembl 数据库ID是ENSG00000188976

 

其余的看头文件自己慢慢理解:

"Functional annotations: 'Allele | Annotation | Annotation_Impact | Gene_Name | Gene_ID | Feature_Type | Feature_ID | Transcript_BioType | Rank | HGVS.c | HGVS.p | cDNA.pos / cDNA.length | CDS.pos / CDS.length | AA.pos / AA.length | Distance | ERRORS / WARNINGS / INFO'

 

软件安装:

选择最新版软件下载:https://sourceforge.net/projects/snpeff/files/

wget https://sourceforge.net/projects/snpeff/files/snpEff_latest_core.zip

因为是java软件,unzip 解压之后就可以直接使用,当然前提是你有java平台。

1

输入数据:

首先下载用来做注释的数据库:java -jar snpEff.jar download GRCh37.75(自己选择需要的版本)

1

软件下载很快,但是数据库下载就需要一定时间啦,去喝杯咖啡吧。

然后软件本身会提供example文件,里面就是一堆各种各样的vcf数据,而且还提供了运行命令,非常简单(examples.sh) ,这些就是我们的输入数据啦!

运行命令:

运行也很简单:java -Xmx4G -jar snpEff.jar -i vcf -o vcf GRCh37.75 example.vcf > example_snpeff.vcf

指定输入输出格式都是vcf,然后指定刚才下载的必备数据库,然后输入输出文件即可!

也可以调用全路径,如果你写在脚本里面的话!

java -Xmx4G -jar path/to/snpEff/snpEff.jar \

-c path/to/snpEff/snpEff.config \

GRCh37.69 \

path/to/example.vcf > example_snpeff.vcf

 

结果解读:

这个非常复杂,对结果理解了多少,就是我们对软件理解了多少。

具体大家看readme吧,注释信息太多了,按需索取:

  1. chromosome_number_variation
  2. exon_loss_variant
  3. frameshift_variant
  4. stop_gained
  5. stop_lost
  6. start_lost
  7. splice_acceptor_variant
  8. splice_donor_variant
  9. rare_amino_acid_variant
  10. missense_variant
  11. inframe_insertion
  12. disruptive_inframe_insertion
  13. inframe_deletion
  14. disruptive_inframe_deletion
  15. 5_prime_UTR_truncation+exon_loss_variant
  16. 3_prime_UTR_truncation+exon_loss
  17. splice_branch_variant
  18. splice_region_variant
  19. splice_branch_variant
  20. stop_retained_variant
  21. initiator_codon_variant
  22. synonymous_variant
  23. initiator_codon_variant+non_canonical_start_codon
  24. stop_retained_variant
  25. coding_sequence_variant
  26. 5_prime_UTR_variant
  27. 3_prime_UTR_variant
  28. 5_prime_UTR_premature_start_codon_gain_variant
  29. upstream_gene_variant
  30. downstream_gene_variant
  31. TF_binding_site_variant
  32. regulatory_region_variant
  33. miRNA
  34. custom
  35. sequence_feature
  36. conserved_intron_variant
  37. intron_variant
  38. intragenic_variant
  39. conserved_intergenic_variant
  40. intergenic_region
  41. coding_sequence_variant
  42. non_coding_exon_variant
  43. nc_transcript_variant
  44. gene_variant
  45. chromosome

http://snpeff.sourceforge.net/VCFannotationformat_v1.0.pdf

 

12

R包精讲第二篇:如何安装旧版本的包?

既然你点进来看,肯定是有需求咯!
一般来说,R语言自带的install.packages函数来安装一个包时,都是默认安装最新版的。
但是有些R包的开发者他会引用其它的一些R包,但是它用的是人家旧版本的功能,但他自己来不及更新或者疏忽了。
而我们又不得不用他的包,这时候就不得不卸载最新版包,转而安装旧版本包。

Continue reading

15

基因组各种版本对应关系

我是受到了SOAPfuse的启发才想到整理各种基因组版本的对应关系,完整版!!!
以后再也不用担心各种基因组版本混乱了,我还特意把所有的下载链接都找到了,可以下载任意版本基因组的基因fasta文件,gtf注释文件等等!!!
首先是NCBI对应UCSC,对应ENSEMBL数据库:
GRCh36 (hg18): ENSEMBL release_52.
GRCh37 (hg19): ENSEMBL release_59/61/64/68/69/75.
GRCh38 (hg38): ENSEMBL  release_76/77/78/80/81/82.
可以看到ENSEMBL的版本特别复杂!!!很容易搞混!
但是UCSC的版本就简单了,就hg18,19,38, 常用的是hg19,但是我推荐大家都转为hg38
看起来NCBI也是很简单,就GRCh36,37,38,但是里面水也很深!
Feb 13 2014 00:00    Directory April_14_2003
Apr 06 2006 00:00    Directory BUILD.33
Apr 06 2006 00:00    Directory BUILD.34.1
Apr 06 2006 00:00    Directory BUILD.34.2
Apr 06 2006 00:00    Directory BUILD.34.3
Apr 06 2006 00:00    Directory BUILD.35.1
Aug 03 2009 00:00    Directory BUILD.36.1
Aug 03 2009 00:00    Directory BUILD.36.2
Sep 04 2012 00:00    Directory BUILD.36.3
Jun 30 2011 00:00    Directory BUILD.37.1
Sep 07 2011 00:00    Directory BUILD.37.2
Dec 12 2012 00:00    Directory BUILD.37.3
可以看到,有37.1,   37.2,  37.3 等等,不过这种版本一般指的是注释在更新,基因组序列一般不会更新!!!
反正你记住hg19基因组大小是3G,压缩后八九百兆即可!!!
如果要下载GTF注释文件,基因组版本尤为重要!!!
对于ensembl:
变幻中间的release就可以拿到所有版本信息:ftp://ftp.ensembl.org/pub/
对于UCSC,那就有点麻烦了:
需要选择一系列参数:
2. Select the following options:
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Genes and Gene Predictions
track: UCSC Genes
table: knownGene
region: Select "genome" for the entire genome.
output format: GTF - gene transfer format
output file: enter a file name to save your results to a file, or leave blank to display results in the browser
3. Click 'get output'.
 现在重点来了,搞清楚版本关系了,就要下载呀!
UCSC里面下载非常方便,只需要根据基因组简称来拼接url即可:
或者用shell脚本指定下载的染色体号:
for i in $(seq 1 22) X Y M;
do echo $i;
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
## 这里也可以用NCBI的:ftp://ftp.ncbi.nih.gov/genomes/M_musculus/ARCHIVE/MGSCv3_Release3/Assembled_Chromosomes/chr前缀
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fasta;
done
rm -fr chr*.fasta
12

bioconductor中文社区招募站长

我已经构建好bioconductor中文社区的雏形,大家可以进去看看!

写在前面

突然发现我的bioconductor.cn这个域名都快要过期了!

哈哈,才想起一年前的计划到现在还没开始实施,实在不像我的风格,可能是水平到了一定程度吧,很多初级工作不像以前那样事无巨细的把关了。正好,借这个机会找几个朋友帮我一起完成这个bioconductor中文社区计划!

 

Continue reading

18

生物信息学分析过程中常见文件格式

刚开始接触生物信息学的时候我也很纠结什么fastq,fastq,sam,bam,vcf,maf,gtf,bed,psl等等,甚至还有过时了的NCBI,ENSEMBL格式,如果是我刚开始 学的时候,我倒是很愿意把他们全部搞透彻,写详细的说明书,但是现在成长了,这些东西感觉很low了,正好我看到了一篇帖子讲数据格式的收集大全,分享给大家,希望初学者能多花点时间好好钻研!

https://www.biostars.org/p/55351/

每种文件格式的定义,都是有它的道理的,大部分是因为一个比较流行的软件,少量的数据格式是因为国际组织广泛认可而流行的
27

gene的symbol与entrez ID并不是绝对的一一对应的

很多时候,我们都无法确定到底是基于symbol来进行分析,还是基于entrez ID,当我们要进行ID转换的时候也想当然的以为它们的一一对应的, 但是最近我写了一个脚本来分析CCLE的数据的时候,发现其实有一些特例:

suppressMessages(library(org.Hs.eg.db))

all_symbol=mappedkeys(org.Hs.egSYMBOL2EG)

all_EGID =mappedkeys(org.Hs.egSYMBOL)
tmp=as.list(org.Hs.egSYMBOL2EG[all_symbol])
#tmp=as.list(org.Hs.egSYMBOL[all_EGID ])
tmp=unlist(lapply(tmp,length))
tmp=tmp[tmp>1]
as.list(org.Hs.egSYMBOL2EG[names(tmp)])
有多个entrez ID对应一个symbol的现象出现,但是没有一个symbol对应多个entrez ID的现象。而且entrez ID也会过期!
$CSNK1E
[1] "1454"      "102800317"
$HBD
[1] "3045"      "100187828"
$RNR1
[1] "4549" "6052"
$RNR2
[1] "4550" "6053"
$SFPQ
[1] "6421"   "654780"
$TEC
[1] "7006"      "100124696"
$MEMO1
[1] "7795"  "51072"
$KIR3DL3
[1] "115653"    "100133046"
$MMD2
[1] "221938"    "100505381"
$`LSAMP-AS1`
[1] "100506708" "101926903"
通过下面的链接可以看到具体情况
14

用broad出品的软件来处理bam文件几次遇到文件头错误

报错如下:ERROR MESSAGE: SAM/BAM file input.marked.bam is malformed: SAM file doesn't have any read groups defined in the header.  The GATK no longer supports SAM files without read groups !

有些人遇到的是bam的染色体顺序不一样,还有可能是染色体的名字不一样,比如>1和>chr1的区别,虽然很傻,但是遇到这样问题的还不少!
还有一些人是遇到基因组没有dict文件,也是用picard处理一下就好。

大部分人是在GATK遇到的,我是在RNA-SeQC遇到的,不过原理都是一样的。
都是因为做alignment的时候并未添加头信息,比如:
bwa samse ref.fa my.sai my.fastq > my.sam
samtools view -bS my.sam > my.bam
samtools sort my.bam my_sorted
java -jar ReordereSam.jar I=/path/my_sorted.bam O=/path/my_reordered.bam R=/path/ref.fa
通过这个代码可以得到排序好的bam,但是接下来用GATK就会报错
java -jar GenomeAnalysisTK.jar -T DepthOfCoverage -R /paht/ref.fa -I /path/aln_reordered.bam
就是因为没有头信息,group相关信息,解决方法有两种:
bwa samse -r @RG\tID:IDa\tSM:SM\tPL:Illumina ref.fa my.sai my.fastq > my.sam
java -jar AddOrReplaceReadGroups I=my.bam O=myGr.bam LB=whatever PL=illumina PU=whatever SM=whatever
一种是比对的时候就加入头信息,这个需要比对工具的支持。
第二种是用picard工具来修改bam,推荐用这个!虽然我其实并不懂这些头文件信息是干嘛的, 但是broad开发的软件就是需要!希望将来去读PHD能系统性的学习一些基础知识!

 

11

关于芯片平台GPL15308和GPL570

它们虽然被GEO数据标记了不同的ID号,但是其实都是一种芯片,都是昂飞公司的U133++2芯片,分析过芯片数据的人肯定不会陌生了

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL15308

事实上,这个平台应该是GPL570,但是被CCLE数据库给稍微变通了一下,就给了一个GPL15308的标签,平台主页也写的很清楚,它的探针ID是伪ID,其实就是entrez gene ID

1

本来这个芯片设计的是五万多个探针,最后只剩下了18926个基因
This array is identical to GPL570 but the data were analyzed with a custom CDF Brainarray Version 15, hgu133plus2hsentrezg.
十二 30

用GSEA来做基因集富集分析

how to use GSEA?
这个有点类似于pathway(GO,KEGG等)的富集分析,区别在于gene set(矫正好的基于文献的数据库)的概念更广泛一点,包括了

how to download GSEA ?

what's the input for the GSEA?

说明书上写的输入数据是:GSEA supported data files are simply tab delimited ASCII text files, which have special file extensions that identify them. For example, expression data usually has the extension *.gct, phenotypes *.cls, gene sets *.gmt, and chip annotations *.chip. Click the More on file formats help button to view detailed descriptions of all the data file formats.
实际上没那么复杂,一个表达矩阵即可!然后做一个分组说明的cls文件即可。
主要是自己看说明书,做出要求的数据格式:http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
表达矩阵我这里下载GSE1009数据集做测试吧!
cls的样本说明文件,就随便搞一搞吧,下面这个是例子:
6 2 1
# good bad
good good good bad bad bad
文件如下,六个样本,根据探针来的表达数据,分组前后各三个一组。
clipboard
现在开始运行GSEA!

start to run the GSEA !

首先载入数据
clipboard
确定无误,就开始运行,运行需要设置一定的参数!
clipboard

what's the output ?

输出的数据非常多,对你选择的gene set数据集里面的每个set都会分析看看是否符合富集的标准,富集就出来一个报告。
点击success就能进入报告主页,里面的链接可以进入任意一个分报告。
最大的特色是提供了大量的数据集:You can browse the MSigDB from the Molecular Signatures Database page of the GSEA web site or the Browse MSigDB page of the GSEA application. To browse the MSigDB from the GSEA application:
 
有些文献是基于GSEA的:

 

十二 22

根据chrY独有区域的覆盖度及测序深度来判断性别

这个也是基于bam文件来的,判断chrY独有区域的覆盖度及测序深度
首先下载chrY独有区域的记录文件,https://www.familytreedna.com/documents/bigy_targets.txt
然后用samtools depth来统计测序深度,samtools  depth $i |grep 'chr[XY]'
depth统计结果文件如下:
mzeng@ubuntu:/home/jmzeng/gender_determination$ head Sample3.depth
chrX    60085    1
chrX    60086    1
chrX    60087    1
chrX    60088    1
chrX    60089    1
chrX    60090    1
chrX    60091    1
chrX    60092    1
chrX    60093    1
chrX    60094    1
然后我随便写了一个脚本来对测序深度文件进行再统计,统计覆盖度及测序深度

[perl]
open FH,"bigy_targets.txt";
while(<FH>){
 chomp;
 @F=split;
 $all+=$F[2]-$F[1]+1;
 foreach ($F[1]..$F[2]){
  $h{$_}=1;
 }
}
close FH;
open FH,$ARGV[0];
while(<FH>){
 chomp;
 @F=split;
 next unless $F[0] eq 'chrY';
 if (exists $h{$F[1]}){
  $pos++;
  $depth+=$F[2];
 }
}
close FH;
$average=$depth/$pos;
$coverage=$pos/$all;
print "$pos\t$average\t$coverage\n" ;

[/perl]

 

 
这样对那三个样本结果如下:
clipboard
可以看到只有sample4,是覆盖率极低的,而且记录到的pos位点也特别少,所以她是女性!
这里测序深度没有意义。
十一 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

 

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

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