R语言:火山图(Volcano Plot)

两个实验组之间的若干个基因的差异倍数比较,可以通过 Volcano Plot 图形化展示。

QQ截图20180425163433

例如 DEseq2 对转录组测序的样品做差异基因分析之后,会得到一个基因表达差异的表格的。案例使用数据的来源是 GitHub:https://github.com/hbc/fabio-splicing/blob/master/new-samples/2016-07-29_fabio-newlines/summary/control-vs-e7107-deseq2.csv

下载上述数据之后,运行下面的命令即可得到文章起始的图。

setwd("C:\\Users\\FengLei\\Desktop\\scatter-plot")
library(ggplot2)
library(ggthemes)
library(Cairo)

data=read.csv("control-vs-e7107-deseq2.csv", head=T)

### 火山图
data$threshold = as.factor(ifelse(data$pvalue =1, 'Up', ifelse(data$pvalue < 0.05 & data$log2FoldChange <= -1, 'Down', 'Not')))
ggplot(data=data, aes(x=log2FoldChange, y=-log10(pvalue), colour=threshold, fill=threshold)) + 
 scale_color_manual(values=c("blue", "grey","red"))+
 geom_point(alpha=0.4, size=1.2) +
# xlim(c(-4, 4)) +
# ylim(c(0, 7.5)) +
 theme_bw(base_size = 12, base_family = "Times") +
 geom_vline(xintercept=c(-1,1),lty=4,col="grey",lwd=0.6)+
 geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)+
 theme(legend.position="right",
 panel.grid=element_blank(),
 legend.title = element_blank(),
 legend.text= element_text(face="bold", color="black",family = "Times", size=8),
 plot.title = element_text(hjust = 0.5),
 axis.text.x = element_text(face="bold", color="black", size=12),
 axis.text.y = element_text(face="bold", color="black", size=12),
 axis.title.x = element_text(face="bold", color="black", size=12),
 axis.title.y = element_text(face="bold",color="black", size=12))+
 labs(x="log2 (Fold Change)",y="-log10 (P-value)",title="Volcano Plot")

注意 R 语言中 ifelse() 函数的使用方法。其基本语法格式如下:
if(con,statement1,statement2)
con 是逻辑条件,当逻辑条件的值为 TRUE 时,则输出 statement1 的值,否则输出 statement2 的值。
也就是说 ifelse() 函数只能输出两个类型的值,然而本文使用的转录组案例中,组间基因变化有三种情况,分别为上调、下调与无明显变化,这里可以使用 ifelse() 的嵌套结构来实现三种判断值的输出。

# ifelse() 函数使用例1
 > x y # 如果向量 x 中的元素值非0,就输出0,否则输出1
 > y
 [1] 0 1 0 1 0 0 1 1
 >
 > # ifelse() 函数使用例2: 嵌套结构的应用
 > x y0, 2*x-1, ifelse(x==0,0,3*x-10))
 > # 如向量 x 中元素值大于0,等于0和小于0,分三种情况做转换
 > y
 [1] 5 19 0 -13 -40
 >

参考资料

[1] ggplot 与火山图:https://www.tanboyu.com/ggplot2-for-volcano.html

[2] R 语言 ifelse() 函数使用:http://www.biye5u.com/article/R/2017/6352.html

Advertisements
发表在 Bioinformatics

mitogenome (unfinished)

 

tRNA

tRNA sequences  整理成 fasta 格式

通过 blastn mapping  到 mitogenome,输出 format=6。

编程序转换成 features table 格式的文件。

注意:

1)部分 tRNA 有多个 copies,例如 trnC(gca) 与 trnfM(cat) 各有两个copies;

2)部分 tRNA 具有同样的 anticodon,但是其名称和序列不同,例如 trnI(cat) 和 trnM(cat);

3)tRNA都没有内含子,所以可以较为简单的转换格式;

4)anticodon 如何自动注释?

 

 

发表在 Bioinformatics

超几何分布与基因功能富集分析

GO功能显著性富集分析给出与基因组背景相比,在差异表达基因中显著富集的GO功能条目,从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向Gene Ontology数据库(http://www.geneontology.org/)的各个term映射,计算每个term的基因数目,然后应用超几何检验,找出与整个基因组背景相比,在差异表达基因中显著富集的GO条目,其计算公式为

hypergeometric

其中,N为所有基因中具有GO注释的基因数目;n为N中差异表达基因的数目;M为所有基因中注释为某特定GO term的基因数目;m为注释为某特定GO term的差异表达基因数目。计算得到的p-value通过Bonferroni校正之后,以corrected p-value≤0.05为阈值,满足此条件的GO term定义为在差异表达基因中显著富集的GO term。通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。

例如,如下的结果:

QQ截图20180208115638

P值计算方法为(在R语言程序中调用phyper命令):

phyper(66-1, 2775, 20269-2775, 249, lower.tail=F)
[1] 5.546399e-08

Q值计算方法:

Q值则需要考虑到全部差异表达基因的P值。引用他人的解释:设总共有m个候选基因,每个基因对应的p值从小到大排列分别是p(1),p(2),…,p(m),则若想控制fdr不能超过q,则只需找到最大的正整数 i,使得 p(i)<= (i*q)/m.然后,挑选对应p(1),p(2),…,p(i)的基因做为差异表达基因,这样就能从统计学上保证fdr不超过q。

例如

> p<-c(0.0003,0.0001,0.02)
> p.adjust(p,method="fdr",length(p))
[1] 0.00045 0.00030 0.02000

 

注:

一、什么是超几何分布

超几何分布的公式为:

p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k)
for x = 0, …, k.
其中, m 是桶里面白球的个数, n 是黑球的个数, k 是从桶中随机取出的球数, x 是取出球
中白球的个数.

Fisher‘s Exact Test就是依据超几何概率分布得到的。

二、计算公式

公式:计算P值的代码也可以自己写:C(k,M)*C(n-k,N-M)/C(n,N)

R语言的代码:

dhyper(x, m, n, k, log = FALSE)
 计算某一个点的p值
phyper(q, m, n, k, lower.tail = TRUE, log.p = FALSE)
 计算一个分布的p值,默认下计算P(X<=x)
qhyper(p, m, n, k, lower.tail = TRUE, log.p = FALSE)
 计算某一个p值对应的分位数
rhyper(nn, m, n, k)
 按超几何分布给出nn的可能的模拟结果值

三、运用案例

(见前文)

四、编程

参考:

http://www.chenlianfu.com/?p=1122

http://www.biotrainee.com/thread-749-1-1.html

http://fhqdddddd.blog.163.com/blog/static/1869915420169212398814

http://blog.sina.com.cn/s/blog_73e6125c01012412.html

 

 

 

发表在 Bioinformatics

安装vcftools

configure的过程报错:Zlib找不到。

实际上已经安装zlib。

解决办法

export ZLIB_LIBS='-L/home/fenglei/local/lib -lz'
export ZLIB_CFLAGS=-I/home/fenglei/local/include
echo $ZLIB_LIBS $ZLIB_CFLAGS
./configure --prefix=/home/fenglei/local/
make
make install

成功安装。

参考资料:https://github.com/vcftools/vcftools/issues/29

发表在 Bioinformatics

[miRNA analysis]安装MiRPlant,patman

安装miRPlant程序,先需要安装patman(一个比对small RNA sequence的程序),但是安装patman的过程中遇到错误了。

[fenglei@localhost PatMaN]$ make install prefix=/home/fenglei/local/
g++ -DVERSION="\"1.2\"" -Wall -O3 -funroll-loops -DNDEBUG -march=k8 -c -o prefix_tree.o prefix_tree.cpp
In file included from /home/fenglei/local/include/assert.h:5:0,
 from /home/fenglei/local/include/c++/6.3.0/cassert:44,
 from prefix_tree.cpp:12:
/home/fenglei/local/include/except.h:15:32: error: conflicting declaration ‘typedef struct Except_Frame_T* Except_Frame_T’
 typedef struct Except_Frame_T *Except_Frame_T;
 ^~~~~~~~~~~~~~
/home/fenglei/local/include/except.h:15:16: note: previous declaration as ‘struct Except_Frame_T’
 typedef struct Except_Frame_T *Except_Frame_T;
 ^~~~~~~~~~~~~~
/home/fenglei/local/include/except.h:17:18: error: field ‘prev’ has incomplete type ‘Except_Frame_T’
 Except_Frame_T prev;
 ^~~~
/home/fenglei/local/include/except.h:16:8: note: definition of ‘struct Except_Frame_T’ is not complete until the closing brace
 struct Except_Frame_T {
 ^~~~~~~~~~~~~~
make: *** [Makefile:25: prefix_tree.o] Error 1

根据错误信息提示,查找下面的文件

/home/fenglei/local/include/except.h

这是GMAP编译过程中生成的文件,于是进入gmap安装目录,执行make uninstall进行卸载。发现except.h文件已经消失。

随后返回patman安装目录,执行make install –prefix=/home/fenglei/local 就顺利安装了。

 

 

发表在 Bioinformatics

Gene Ontology analysis for DEGs of Arabidopsis

try http:// if https:// URLs are not supported
 source("https://bioconductor.org/biocLite.R")
 biocLite("clusterProfiler")
 biocLite("DOSE")
 biocLite("tibble")
 library(clusterProfiler)
 biocLite("topGO")
 library(topGO)
 biocLite("org.At.tair.db")
 library(org.At.tair.db)

a=read.table("diff.gene.table", head=T, sep="\t")a=read.table("diff.gene.table", head=T, sep="\t")b=a[,1]
keytypes(org.At.tair.db)  ## 看该数据库支持哪些基因名称类型,例如拟南芥支持AT1G01110 就是keytype="TAIR"
ego <- enrichGO(gene          = b, keyType = "TAIR",                OrgDb         = org.At.tair.db,                ont           = "CC",                pAdjustMethod = "BH",                pvalueCutoff  = 0.05,                qvalueCutoff  = 0.05, readable      = TRUE)
barplot(ggo, drop=TRUE, showCategory=12)dotplot(ego)
write.table(as.data.frame( ego@result), file="test_CC.txt")
# KEGG over-representation testkk <- enrichKEGG(gene         = b,                 organism     = 'ath',                 pvalueCutoff = 0.05)
 write.table(as.data.frame(kk@result), file="test_kk.txt")

# KEGG Gene Set Enrichment Analysiskk2 <- gseKEGG(geneList     = b,               organism     = 'ath',               nPerm        = 1000,               minGSSize    = 120,               pvalueCutoff = 0.05,               verbose      = FALSE)head(kk2)
mkk <- enrichMKEGG(gene = b,                   organism = 'ath')

 

发表在 Bioinformatics

根据参考基因组对contigs进行排列

Multiple programs have been developed for reference-assisted chromosome assembly: Bambus [10], BACCardI [11], Projector2 [12], OSLay [13], ABACAS [14], MeDuSa [15], AlignGraph [16], Ragout [17], SyMap [18] and RACA [19]. Most of the listed tools were designed for bacterial or small genomes. For example, ABACAS is a convenient bacterial genome contiguation tool that may also be used for small eukaryotic genomes such as Saccharomyces cerevisiae (12.1 mega base pairs). However, ABACAS is not efficiently scaled to use with the large genomes typical of vertebrate species.

ABACAS对小型基因组很实用。但是注意:contig连接处是由99个N连接起来的,实际上contig的末端室友重叠关系的。

针对大型基因组,可以使用Chromosomer: a reference-based genome arrangement tool for producing draft chromosome sequences

发表在 Bioinformatics