差异分析
一、转录本表达矩阵构建
1.转录本表达定量
align_and_estimate_abundance.pl \ #进行转录本表达定量
--transcripts Trinity.fasta \#输入文件是Trinity组装出的Trinity.fasta用作参考基因组
--seqType fq \ #读取的原始数据的数据类型是fastq
--left fa_1.fastq.gz --right fa_2.fastq.gz \#质控后的clean.data比对定量
--est_method RSEM \#进行表达定量的计算软件RSEM
--aln_method bowtie2 \RSEM调用比对工具为bowtie2
--trinity_mode \
--prep_reference \
--output_dir rsem \
--thread_count 70
调用脚本得到的数据
bowtie2.bam #bowtie2 生成的 bam文件
bowtie2.bam.for_rsem.bam #用于RSEM计算的 bam文件
bowtie2.bam.ok
RSEM.genes.results #基于基因的RSEM Reads Count
RSEM.isoforms.results #基于转录本的 RSEM Reads Count
RSEM.isoforms.results.ok
RSEM.stat
├── RSEM.cnt
├── RSEM.model
└── RSEM.theta
由于我们在每一个文件夹中的Reads count 没有经样本间的均一化，因此需要做一个样本均一化，构建转录本-基因表达矩阵并得到不同样本中的均一化表达数据TMM
2.将不同样本比对结果整理成一个文件并进行均一化处理
find * -name '*.isoforms.results'> quant.file
#
采用了find命令将子文件夹中的RSEM.isoforms.results
基因表达量结果全部查找出来并将路径放到一个文件中（quant.file）用以后期使用
RSEMmatrix.sh abundance_estimates_to_matrix.pl 整理文件用于差异分析
abundance_estimates_to_matrix.pl \
--est_method RSEM \ #由于是对RSEM的结果进行矩阵构建，因此需要说明这个
--gene_trans_map \
--name_sample_by_basedir \
--quant_files quant.file
代码注释
#--est_method <string> RSEM|eXpress|kallisto|salmon (用什么软件分析就填什么软件)
#--gene_trans_map <string> 基因与转录本对应文件(如果不知道填none)
#可以从gtf中提取转录本信息，最后做成的文件要转录本-基因一一对应，转录本一定要在前边，两者用\t分隔。
#--cross_sample_norm <string>TMM|UpperQuartile|none (default: TMM)
#--name_sample_by_basedir 一定要填不然会报错
#--basedir_index <int> default(-2)
#--out_prefix <string> 前缀随意命名，默认为RSEM
#--quant_files <string> 可以是文件夹或者手动输入需要建立矩阵的文件.
样本均一化RSEMmatrix.sh abundance_estimates_to_matrix.pl输出的内容
RSEM.isoform.TPM.not_cross_norm.TMM_info.txt
RSEM.gene.TPM.not_cross_norm.TMM_info.txt
RSEM.gene.TPM.not_cross_norm.runTMM.R
RSEM.isoform.TPM.not_cross_norm.runTMM.R
RSEM.gene.TPM.not_cross_norm
RSEM.gene.counts.matrix
RSEM.gene.TMM.EXPR.matrix
RSEM.isoform.TPM.not_cross_norm
RSEM.isoform.counts.matrix
RSEM.isoform.TMM.EXPR.matrix
主要用以下两个文件进行分析
- '.counts.matrix' 文件用于后期的差异表达分析
- '.TMM.EXPR.matrix'文件可以用于其他基因表达的分析
3.counts.matrix 后续差异性统计分析
1.建立生物性重复对应文件
B25 B251
B25 B252
R25 R251
R25 R252
W25 W251
W25 W252
需要注意的是：samples.txt中的名字需要和matrix中的名字一致，否则没办法识别
2.gene.matrix.count文件 或者isoform.counts.matrix文件
用edgeR进行差异性分析
用R进行分析
安装软件
source("http://bioconductor.org/biocLite.R")
bioLite('limma')
bioLite('edgeR')
bioLite('DESeq2')
bioLite('ctc')
bioLite('Biobase')
bioLite('gplots')
bioLite('ape')
有生物学重复代码
run_DE_analysis.pl \ #调用脚本
--matrix isoforms.counts.matrix \ #abundance_estimates_to_matrix.pl生成的矩阵文件
--method edgeR \ #调用的分析软件
--samples_file samples.isoforms.txt \ #生物学重复信息文件
--output DEEdgeRout #结果输出文件夹名
生成结果文件
Rscript 分析脚本
cout_matrix 两两比对数据
DE_results 比对结果
PDF 可视化图片
4.差异基因注释
1.绘制差异基因热图和样本分布
调用analyze_diff_expr.pl
analyze_diff_expr.pl \
--matrix isoforms.TMM.EXPR.matrix \
--samoles.isoforms.txt \
-P 0.001 \
-C 2
结果文件
1. 两两比较的样本信息（.samples）
2. 样本比较后的基因表达情况上调还是下调的文件以及整合在一起的基因（.subset）
3. 差异基因在不同样本中分布数量（DE_feature_counts.P0.001_C2.matrix）
4. 差异基因表达量矩阵（diffExpr.P0.001_C2.matrix）
5. 差异基因对数转换后的表达矩阵 （diffExpr.P0.001_C2.matrix.log2.centered.dat）
6. 出的两个图（PDF）
2.差异基因筛选
通过edgeR这个R包对之前得到的表达量矩阵进行了差异表达验证，其中两个参数很重要一个是FDR和FoldChange，通过设定这两个值可以筛选出想要的差异基因
## 安装limma包
source("http://bioconductor.org/biocLite.R")
biocLite("limma")
## 读入limma包
library(limma)
#读取同一目录下的所有文件
path<-"D:/expr"#文件目录
str_path="D:/txt"
fileNames<-dir(path) #获取该路径下的文件
filePath<-sapply(fileNames, function(x){
paste(path,x,sep='/')}) #生成读取文件路径
for(i in 1:length(fileNames)){
filepath = filePath[i]
strpath = paste(str_path,sub("csv","txt",fileNames[i]),sep='/')
data = read.csv(filepath, header=T, row.names=1)
tag = read.table(strpath, header=F, stringsAsFactors=F)#读取数据，结果为list
summary(tag)##get the file name one by one
##read the expression data and sample tag
out_name = sub(".csv","_DEG.txt",fileNames[i])
# group names for all samples
sml <- c()
gsms<-tag
for (i in 1:nchar(gsms)) { sml[i] <- substr(gsms,i,i) }
sml <- paste("G", sml, sep="") # set group names
design <- model.matrix(~0+factor(sml))
colnames(design) <- c("GA","GB") #GA-case;GB-control
fit <- lmFit(data, design)
cont.matrix <- makeContrasts(GA-GB, levels=design)#
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
tT <- topTable(fit2, adjust="fdr", sort.by="B", number=length(fit2$p.value))
tT_filter <- tT[tT$adj.P.Val<=0.05,] # 以小于0.05筛选出p值显著的基因
tT_filter <- tT_filter[abs(tT_filter$logFC)>=1,] # 以logFC>=1筛选出显著差异表达的基因
#output the DEG data
write.table(tT_filter, out_name, row.names=T, sep="\t")
}
富集分析
library(dplyr)
library("clusterProfiler") #富集通路分析
library("org.Hs.eg.db") #基因名称转换
require(DOSE)
library(enrichplot) #富集分析
library(GSEABase) #富集分析
#GENE_list= row.names(tT_filter) #取基因列表
#result = data.frame(tT_filter,GENE_list)
#output the DEG data
#将基因由ENTREZID/SYMBOL转换成ENSEMBL
ALLgenes = keys(org.Hs.eg.db,keytype = "ENSEMBL")
list = select(org.Hs.eg.db,keys=ALLgenes,columns = c('ENTREZID',"SYMBOL"), keytype="ENSEMBL")
go_gene <- rownames(tT_filter)
go_gene_list = list[match(go_gene,list[,'ENSEMBL']),]
result = data.frame(tT_filter,go_gene_list)
result.sort <- arrange(result, desc(logFC))
write.table(result.sort, out_name, row.names=T,sep="\t")#
entrezIDs = go_gene_list[,'ENTREZID']
# GO 通路富集1
ego = enrichGO(OrgDb = 'org.Hs.eg.db',gene = entrezIDs,ont = "BP",pvalueCutoff = 0.05, readable= TRUE)
write.csv(summary(ego),sub(".txt","_GO.csv",out_name),row.names =F)
dsamp <- diamonds[sample(nrow(diamonds), 1000), ]
pdf(sub(".txt","_goplot.pdf",out_name))
goplot(ego)
dev.off()
pdf(sub(".txt","_barplot20.pdf",out_name))
barplot(ego,showCategory = 20)
dev.off()
# GO 通路富集2
enrichGO <- enrichGO(gene = entrezIDs, universe=universe, organism = "human", ont = "CC", pvalueCutoff = 0.05, readable= TRUE)
#barplot(enrichGO,drop=TRUE,showCategory = 1000)
write.csv(summary(enrichGO),sub(".txt","_GO.csv",out_name),row.names =F)
# kEGG 通路富集1
ekk <- enrichKEGG(gene = entrezIDs,organism = 'hsa',pvalueCutoff = 0.05)
write.csv(summary(ekk),sub(".txt","_KEGG.csv",out_name),row.names =F)
pdf(sub(".txt","_dotplot_KEGG_30.pdf",out_name),width = 10)
dotplot(ekk,showCategory = 30)
dev.off()
# kEGG 通路富集2
library("curl")
ekk <- enrichKEGG(gene=entrezIDs, universe=universe, organism = "human", pvalueCutoff=0.05, readable= TRUE)
if(!is.na(ekk)){
write.csv(summary(ekk),sub(".txt","_KEGG.csv",out_name),row.names =F)
gene1 = result.sort$logFC
names(gene1) = result.sort$ENTREZID#SYMBOL
gene1 = sort(gene1,decreasing = T)
# gene_concept network
cnet_outfile1 = sub(".txt","_cnetplot_GO.pdf",out_name)
cnet_outfile1 = paste('/Users/xinmiaoyan/Desktop/Elva/paper/Drug resistance/GDSC/analysis/IC50_data/expression_data/DEG_GeneSymbol_result/GO_CnetPlot',cnet_outfile1,sep='/')
pdf(sub(".txt","_cnetplot_GO.pdf",out_name))
cnetplot(ego, foldChange = gene1)
dev.off()
PDF(cnet_outfile1)
# GSEA富集分析
kk = gseKEGG(gene1,nPerm = 1000,pvalueCutoff = 0.5)
pdf(sub(".txt","_ridgeplot_p5.pdf",out_name))
ridgeplot(kk)
dev.off()
pdf(sub(".txt","_gseaplot_p5.pdf",out_name))
gseaplot(kk,geneSetID = 1,title = kk$Description[1])
dev.off()
'''
'''
d='./MsigDB/symbols'
gmts <- list.files(d,pattern = 'all')
gmts
'''
'''
gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
c5 <- read.gmt(gmtfile)#前面俩行大概就是来选择GSEA的数据库用的，c5.cc.v5.0.entrez.gmt这里他事例用的是GO的cc数据库，这个可以调整的，参考GSEA的网站来选择
egmt <- enricher(gene1, TERM2GENE=c5)
write.table()write.csv(summary(egmt),"Gsea-enrich.csv",row.names =F)
