17

一文学会WGCNA分析

WGCNA 分析

基本概念

WGCNA其译为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块(module),并探索基因网络与关注的表型之间的关联关系,以及网络中的核心基因。

适用于复杂的数据模式,推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有:不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。 Continue reading

12

perl模块安装大全

今天又有小伙伴微信问我perl模块安装的问题,因为ENSEMBL发布的大多数数据库以及软件都是基于perl的,尤其是分量很重的VEP,所以即使你再如何如何的讨厌perl,也不得不与之打交道。

这种细节问题问我,我当然无法直接给出答案咯。毕竟,我的知识积累都不是靠死记硬背的。所以需要取回过头查看一下我的博客,才意识到,我竟然已经写了7篇教程,关于perl的模块。目录如下:

首先需要自己确定已经安装了哪些模块,都安装在哪里?还有新的模块需要安装到哪里?
然后再学习如何安装新的模块。

Continue reading

05

文献阅读笔记-2014-K27M-H3.3-DIPG

DIPG和 adult glioblastomas (GBMs) 很大区别,不应该用同样的治疗方式。

测序策略是:

We integrated deep sequencing analysis of 36 tumor-normal pairs (20 whole-genome sequencing (Illumina HiSeq 2000) and 16 whole-exome sequencing (Applied Biosystems SOLiD 5500xl)) with comprehensive methylation (28 DIPGs; Illumina Infinium450k methylation array), copy number (45 DIPGs; Affymetrix SNP6.0) and expression (35 DIPGs; Illumina HT-12 v4) data (Supplementary Table 1).  Continue reading

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