十二 15

制作自己的gene set文件给gsea软件

熟悉GSEA软件的都知道,它只需要GCT,CLS和GMT文件,其中GMT文件,GSEA的作者已经给出了一大堆!就是记录broad的Molecular Signatures Database (MSigDB) 已经收到了18026个geneset,但是我奇怪的是里面竟然没有包括cancer testis的gene set,MSigDB的确是多,但未必全,其实里面还有很多重复。而且有不少几乎没有意义的gene set。那我想做自己的gene set来用gsea软件做分析,就需要自己制造gmt格式的数据。因为即使下载了MSigDB的gene set,本质上就是gmt格式的数据而已:http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats#GMT:_Gene_Matrix_Transposed_file_format_.28.2A.gmt.29 Continue reading

十二 15

CpG Islands记录文件下载的4种方式

这个也是读者来信最多的,关于基因组某些区域的起始终止坐标的下载问题,genomic feature的问题,一般是gtf文件或者bed文件,比如人类hg19上面的所有外显子的坐标记录文件,所有基因的坐标记录文件,所有lncRNA,rRNA等等,我这里拿CpG Islands记录文件下载的4种方式举例子给大家说明一下: Continue reading

十二 11

gene symbol 中的奇怪开头基因

这本是我为论坛的基础板块写的一个基础知识点,但是浏览量实在有限,不忍它蒙尘,特在博客重新发布一次!原帖见:http://www.biotrainee.com/thread-511-1-1.html

gene symbol 是非常官方的,由HUGO 组织负责维护,有专门的数据库HGNC database of human gene names | HUGO
以前分析数据的时候,有一些基因的symbol很奇怪,让我百思不得其解,比如
C orf 系列基因,
HS.系列基因,
KRTAP系列基因,
LOC系列基因,
MIR系列基因,
LINC系列基因
它们往往一个系列,就有好几百个基因;
C12orf44; Chromosome 12 Open Reading Frame 44;  这个是C orf系列基因的意思
MIR系列基因应该是 miRNA相关的基因
LINC系列基因应该就是long intergenic non-protein coding RNA
LOC系列基因,是非正式的,推定的,日后可能被更合适的名字替代
我这里做好了所有的基因对应关系,去生信菜鸟团QQ群里下载吧,共47938个基因的symbol和entrez gene id还有name,还有alias的对应!

1
还有一些RNA基因,根本就没有symbol,比如:CTA/B/C/D系列的
Aliases for ENSG00000271971 Gene
Quality Score for this RNA gene is 1
Aliases for ENSG00000271971 Gene
CTD-2006H14.2 5
External Ids for ENSG00000271971 Gene
Ensembl: ENSG00000271971
还有,如果你看到HS.开头的基因,它是unigene的ID了,已经不再是symbol啦。

十二 11

用R获取芯片探针与基因的对应关系三部曲-NCBI下载对应关系

这是系列文章,请先看:

用R获取芯片探针与基因的对应关系三部曲-bioconductor

ncbi现有的GPL已经过万了,但是bioconductor的芯片注释包不到一千,虽然bioconductor可以解决我们大部分的需要,比如affymetrix的95,133系列,深圳1.0st系列,HTA2.0系列,但是如果碰到比较生僻的芯片,bioconductor也不会刻意为之制作一个bioconductor的包,这时候就需要自行下载NCBI的GPL信息了,也可以通过R来解决:

##本质上是下载一个文件,读进R里面,然后解析行列式,得到芯片探针与基因的对应关系,看下面的代码,你就能理解了。 Continue reading

十二 01

GSEA的统计学原理试讲

GSEA这个java软件使用非常方便,只需要根据要求做好GCT/CLS格式的input文件就好了。我以前也写个用法教程:

但说到统计学原理,就有点麻烦了,我试着用自己的思路阐释一下:
假设芯片或者其它测量方法测到了2万个基因,那么这两万个基因在case和control组的差异度量(六种差异度量,默认是signal 2 noise,GSEA官网有提供公式,也可以选择大家熟悉的foldchange)肯定不一样,那么根据它们的差异度量,就可以对它们进行排序,并且Z-score标准化,在下图的最底端展示的就是

Continue reading

十二 01

吐血推荐snpedia数据库,非常丰富的snp信息记录

正好,我拿到了自己的全基因组测序数据,而前些天看到朋友圈推送的文章提到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因,我就想顺便看看自己在这两个位点的变异情况。一般的流程是先找完变异位点,然后用vep/snpEFF对变异位点进行注释,然后看看有没有这两个位点。但我仅仅是想查看这两个位点,所以我会根据它的rsID来找到它的基因组坐标,再直接call这个位置的变异情况。以前我都是用dnSNP来查看rsID的基因组坐标的,
mkdir -p ~/annotation/variation/human/dbSNP
cd ~/annotation/variation/human/dbSNP
## https://www.ncbi.nlm.nih.gov/projects/SNP/
## ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh38p2/
## ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/
nohup wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz &
wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz.tbi

Continue reading

十一 25

用BioNet这个bioconductor包来找 maximal-scoring subgraph

## 此包是为了解决一个难题: maximal-scoring subgraph (MSS) problem ,在一个巨大的复杂网络里面找到significantly differentially expressed subnetworks,就是说,得到了几百个差异基因,去PPI数据库做网络图的时候,发现还是巨大无比,所以需要用这个包来精简我们的网络图。
heuristically的中文意思:启发性地
## 而这个R包可以整合多种数据结果来给一个网络打分,
它整合了PPI网络分析和寻找功能模块的需求。
重点就是根据一个"igraph" or "graphNEL"对象和打分来找最大的MSS
subnet <- subNetwork(dataLym$label, interactome)
module <- runFastHeinz(subnet, scores)
plotModule(module, scores=scores, diff.expr=logFC) #这个就是精简后的我们的网络图。
其实另外一个函数也有类似的功能,dNetFind https://rdrr.io/cran/dnet/man/dNetFind.html

Continue reading

十一 24

cytoscape五步曲之二:在cytoscape里面生成网络图

通过上一讲大家应该明白了,网络图是为了展现分子之间的连接关系的,并不是一定要用cytoscape来做,只需要根据连接关系给我们的所有点安排一个坐标,然后把相应的线连接起来即可!那么既然我们要学习cytoscape,肯定是要用cytoscape做好第一步,就是根据输入数据来做网络图。
可以先了解一下cytoscape定义好的输入数据,
http://wiki.cytoscape.org/Cytoscape_User_Manual/Network_Formats 当然,其实木有意义!因为我们不可能拿到cytoscape的输入文件(cys格式的),除非是你朋友传给你的。我们肯定是根据txt.csv等分割的文本文件来做网络图。

Continue reading

十一 24

cytoscape五步曲之一:明白什么是网络图

想了想还是写一个系列教程吧,问的朋友也太多了,主要是因为cytoscape跟python一样,经历了从2到3的进化阵痛过程,而且进化的面目全非了!!!很多人拿着2.x的说明书教程,视频,然后下载的却是3.x版本的cytoscape,真可怕!!!
已经从两万个芯片探测到的基因里面找到了近千个差异基因了,对它们做了GO/KEGG分析还是抓不住重点,看到文献说可以用PPI数据库做network analysis之后找hub基因,也也许可以说明一些问题!
提到 network analysis ,我想起来我以前总结过 R语言画网络图的三部曲,里面讲到过网络分析的基本原理!

Continue reading

十一 24

cytoscape五步曲之三:安装各种插件

软件安装我就不多说了,直接去官网下载即可,请务必下载3.x版本,我讲的是 最新版教程!

本次讲解如何给cytoscape安装插件,cytoscape本身是一个平台,学者可以在上面开发各种各样功能的插件实现不同的分析需求,类似于R语言这个平台,人们在上面安装包一样。R里面如何安装包我博客讲了4次,基本上看完的人都会懂。而cytoscape不一样,它的插件安装非常简单!非常简单!非常简单!

你只需要去cytoscape的APP中心找到包,如果你打开了cytoscape的界面,那么网页就会有install的字样,非常显眼,点击就自动安装了,这个时候会安装到

C:\Users\jimmy1314\CytoscapeConfiguration\3\apps\installed 这个目录!!~ 在你的电脑里面 jimmy1314 不一样

如果你这个时候并没有打开cytoscape的界面,那么网页就会有download的字样,也是非常显眼,点击就可以下载, 下载之后你需要自己把下载的jar文件放到cytoscape的安装路径,一般默认是

C:\Program Files\Cytoscape_v3.3.0\apps

最后,cytoscape提供了APP中心,就跟苹果手机安卓手机安卓软件一样,直接在cytoscape软件的菜单栏app中心就可以点击安装!

我要说的就是这么多了,我安装了十几个插件了,都没有什么问题,如果大家有遇到安装不了的,随时报告我,我来更新教程!联系jmzeng1314@163.com 

下面的链接选择性观看:

http://wiki.cytoscape.org/Cytoscape_3/UserManual

十一 23

quantile normalization到底对数据做了什么?

提到normalization很多人都烦了,几十种方法,而对于芯片或者其它表达数据来说,最常见的莫过于quantile normalization啦。那么它到底对我们的表达数据做了什么呢?首先要么要清楚一个概念,表达矩阵的每一列都是一个样本,每一行都是一个基因或者探针,值就是表达量咯。quantile normalization 就是对每列单独进行排序,排好序的矩阵求平均值,得到平均值向量,然后根据原矩阵的排序情况替换对应的平均值,所以normalization之后的值只有平均值了。具体看下面的图: Continue reading

十一 23

用R的bioconductor里面的stringDB包来做PPI分析

PPI本质上是根据一系列感兴趣的蛋白质或者基因(可以是几百个甚至上千个)来去PPI数据库里面找到跟这系列蛋白质或者基因的相互作用关系!

本次的主角是stringDB,顾名思义用得是大名鼎鼎的string数据库,
本来还以为需要自己上传自己的基因给这个数据库去做分析,没想到他们也开发了R包,主页见: http://www.bioconductor.org/packages/release/bioc/html/STRINGdb.html 而我比较喜欢用编程来解决问题,所以就学了一下这个包,非常好用!
它只需要一个3列的data.frame,分别是logFC,p.value,gene ID,就是标准的差异分析的结果。
然后用string_db$map函数给它加上一列是 string 数据库的蛋白ID,然后用string_db$add_diff_exp_color函数给它加上一列是color。
用string_db$plot_network函数画网络图,只需要 string 数据库的蛋白ID,如果需要给蛋白标记不同的颜色,需要用string_db$post_payload来把color对应到每个蛋白,然后再画网络图。
也可以直接用get_interactions函数得到所有的PPI数据,然后写入到本地,再导入到cytoscape进行画图

Continue reading

十一 23

java环境变量的问题

有篇文章提到了cytoscape,想着一直没用过这个神器对不起我生信大神的称号呀,就下载了准备安装,居然报错了,简直不可思议,因为一直以为它是java软件,一般不需要安装,结果是exe的,只是依赖于java,报错是EXE4J_JAVA_HOME, No JVM could be found on your system,这是个很常见的错误,我 简单搜索了解决方案https://wincrunch.com/exe4j-java-home-no-jvm-could-be-found-on-your-system/ 居然无效,但是里面有句话引起了我的注意,通常64位的window电脑的java是安装在Program Files 而不是Program Files (x86),这才是问题所在,我当初图简单,直接用了JDK来安装JRE,所以导致软件安装目录错误。有非常多的生物信息学软件都依赖与java,比如IGV,GSEA,cytoscape,一般来说window电脑安装好了java之后这些软件都挺好用的。那么关于java问题,我整理了3个:

Continue reading

十一 14

TPM值就是RPKM的百分比嘛!

很久以前就有人问过这个问题啦,虽然目前主流还是用RPKM/FPKM来形容一个基因的表达量。但是既然大家都说TPM更好,我也来探究一下吧!

我不喜欢看公式,直接说事情,我有一个基因A,它在这个样本的转录组数据中被测序而且mapping到基因组了 5000个的reads,而这个基因A长度是10K,我们总测序文库是50M,所以这个基因A的RPKM值是 5000除以10,再除以50,为10. 就是把基因的reads数量根据基因长度和样本测序文库来normalization 。 Continue reading

十一 14

仅仅对感兴趣的基因call variation

有这个需求,是因为我们经常对某些细胞系进行有针对性的设计变异,比如BAF155的R1064K呀,H3F3A的K27呀,那我我们拿到高通量测序数据的时候,就肯定希望可以快速的看看这个基因是否被突变成功了。现在比对几乎不耗费什么时间了,但是得到的sam要sort的时候还是蛮耗费时间的。假设,我们已经得到了所有样本的sort好的bam文件,想看看自己设计的基因突变是否成功了,可以有针对性的只call 某个基因的突变!

1

Continue reading

十一 12

仔细探究picard的MarkDuplicates 是如何行使去除PCR重复reads功能的

本帖紧跟前面的仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的

同样的我们也是分单端和双端测序来看结果,并且比较两个工具的区别!

首先对于那个单端数据,samtools给出的结果是:[bam_rmdupse_core] 25 / 53 = 0.4717 in library Continue reading

十一 12

仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的

在做这个去除PCR重复reads时候必须要明白为什么要做这个呢?WGS?WES?RNA-SEQ?CHIP-SEQ?都需要吗?随机打断测序才需要?特异性捕获不需要?
搞明白了,我们就开始做,首先拿一个小的单端测序数据比对结果来做测试!
samtools rmdup -s tmp.sorted.bam tmp.rmdup.bam
[bam_rmdupse_core] 25 / 53 = 0.4717 in library
我们的测试数据里面有53条records根据软件算出了25条reads都是PCR的duplicate,所以去除了!

Continue reading

十一 07

mysql的table居然有最大列限制

想着把TCGA的RPKM值矩阵表格写入到mysql,然后做一个查询网页给生物学家,我下载的是所有TCGA收集的mRNA表达数据集数据集-GSE62944 ,共9264个癌症样本,和741个正常组织的表达数据。当我想写入癌症表达矩阵的时候,报错了:
Error in .local(conn, statement, ...) :
could not run statement: Too many columns
简单搜索了一下,发现是mysql有最大列的限制,但是我不是很懂计算机,所以没太看明白该如何调整参数使得mysql列限制扩充:http://dev.mysql.com/doc/refman/5.7/en/column-count-limit.html 所以就把癌症表达矩阵根据癌症拆分了,癌种数量如下:

Continue reading