本期使用R的edgeR包对去除了有注释样本的miRNA表达谱进行差异表达分析.

之前已经做了以下几期:

我们已经通过TCGA的API下载了肝癌组织373个cases的425个miRNA表达谱文件和临床数据文件,整合了表达谱文件和临床数据,并去除了有annotations.txt文件的样本.

本期使用R的edgeR包寻找 肝癌 VS 正常 的差异表达miRNA:

1.读入表达谱和样本信息文件:

rm(list=ls())
library(data.table)
library(edgeR)
library(stringi)
# 表达谱文件路劲:
miRNA_exp_path="./TCGA_liver_miRNA/exp_matrix_reads.tsv"
# 样本信息+临床信息文件路径:
cases_annotation_path="./TCGA_liver_miRNA/cases_files_merged.tsv"
# 使用data.table的fread读入表达谱文件(reads count)
miRNA_exp<-fread(miRNA_exp_path,showProgress = T,header = T,sep = "\t")
# 去掉最后一列空行,如果是read.table则不会有这个问题
names(miRNA_exp)[ncol(miRNA_exp)]<-"final_col"miRNA_exp<-miRNA_exp[,!c("final_col")]
# 使用data.table的fread读入样本说明信息
cases_annot<-fread(cases_annotation_path,showProgress = T,header = T)

2.查看样本信息:

# 获得样本信息:
grp<-as.factor(merge(data.table(file_id=names(miRNA_exp)),cases_annot[,c("file_id","case_id","cases_0_samples_0_sample_type")],by="file_id")[,cases_0_samples_0_sample_type])
# 查看样本信息:
table(grp)

之前我们已经去掉了有annotations.txt样本的23个表达谱,现在剩下402样本.结果发现这402个样本中有3个样本来自复发的肝癌,为了避免可能存在的影响,把这3个也去掉吧,反正样本多:

# 发现有3个复发的,直接在表达量矩阵中去掉:
recurrent_tumor_id<-merge(data.table(file_id=names(miRNA_exp)),cases_annot[,c("file_id","case_id","cases_0_samples_0_sample_type")],by="file_id")[,file_id][which(grp=="Recurrent Tumor")] miRNA_exp<-miRNA_exp[,-recurrent_tumor_id,with=F]

再次查看样本信息:

# 再次获得样本信息:
grp<-as.factor(merge(data.table(file_id=names(miRNA_exp)),cases_annot[,c("file_id","case_id","cases_0_samples_0_sample_type")],by="file_id")[,cases_0_samples_0_sample_type])
# 查看样本信息:
table(grp)

3.重新命名表达谱的行名和列名并转为data.frame,开始的时候是使用data.table的fread读入的,因为这个读入和运行速度快,建议大家也慢慢从data.frame转入data.table.

# 保存miRNA名称为列名miRNA_names<-miRNA_exp[,miRNA]
miRNA_exp<-miRNA_exp[,!c("miRNA")]
# 将data.table转为data.frame,因为data.table没有行名,前面使用data.table是因为速度快setDF(miRNA_exp)
# 添加列名
rownames(miRNA_exp)<-miRNA_names
# 拷贝个新的
miRNA_exp_2<-copy(miRNA_exp)
# 使用grp+列数修改列名
grp_chr<-as.character(grp)
grp_chr[which(grp_chr=="Primary Tumor")]<-"T"grp_chr[which(grp_chr=="Solid Tissue Normal")]<-"N"grp<-as.factor(grp_chr)
colnames(miRNA_exp_2)<-paste(grp,1:length(grp),sep="")

4.查看样本之间的关系:

# 查看任意样本表达量直接的关系
pairs(log2(1+miRNA_exp_2[,1:7]), pch=".",lower.panel=NULL)

这个没法看全部399个样本之间的表达量关系,因为太多了画不出来.上面只查看了最前面7个样本的.从上图可以看出样本间的表达量基本都是呈线性的.

5.建立DGEList对象并对数据进行标准化

# 建立一个DGEList对象.
d <- DGEList(counts=miRNA_exp_2, group=grp)
# 使用TMM方法归一化:d <- calcNormFactors(d,method = "TMM")
head(d$samples)

6.样本质量控制:同一组的生物重复应该是聚在一起的,这里即Tumor的样本要聚在一起,Normal的要聚在一起,做法如下:

cols <- as.numeric(d$samples$group) plotMDS(d,col=cols)

从图中可以看出,左下角的N139和N77,中间的N253,N325和N243正常样本并没有聚在一起,而且Tumor样本也有很多比较离散. 并且在右边正常组织和肿瘤组织聚集在了一起.我们先不管,接着往下做.

7.构建模型。

mm <- model.matrix(~-1+grp) head(mm)

8.依据模型估算离散度（dispersion).离散度为生物变异系数（biological coefficient of variation）的平方,衡量样本之间的差异程度,同一组样本之间离散度越小越好.

d <- estimateGLMCommonDisp(d,mm) d <- estimateGLMTrendedDisp(d,mm) d <- estimateGLMTagwiseDisp(d,mm)

使用plotBCV来察看同一组内不同样品的不同表达水平值的方差分析。

对于plotBCV给出的红线，它在0.2~0.4之间都是容易接受的。但如果超过0.5，则需要考虑实验样品之间的差异是否过大了。我们这里得到的变异系数达到了很高的1.28,所以TCGA这399个样本之间的差异太大了,我们先不管,继续往下做.

plotMeanVar(d, show.raw=TRUE, show.tagwise=TRUE, show.binned=TRUE)

从上图可以看出基因水平的差异（variance）随着表达值的升高而升高

plotSmear(d, pair=c("N","T"), ylim=c(-5,5))

从上图可以看出平均表达量和差异表达倍数之间的关系,基本关于Y轴对称,并且可以看出无论表达量高低,绝大部分miRNA是没有差异表达的.

9.接下来先做广义线性拟合（Generalized linear models），而后构建比较结构（construct contrast），最后依据比较结构进行比较，给出差异表达的基因。

f <- glmFit(d,mm) con <- makeContrasts("T-N"=grpT-grpN,levels=colnames(mm)) conlrt <- glmLRT(f,contrast=con) topTags(lrt,20)

这样就得到了每个miRNA的logFC和Pvalue, 将结果保存到文件:

miRNA_DE_result<-topTags(lrt,Inf)$table #获取列名: miRNA_DE_result<-cbind(miRNAs=rownames(miRNA_DE_result),miRNA_DE_result) write.table(miRNA_DE_result,file="./TCGA_liver_miRNA/edgeR_miRNA_DE.tsv",quote = F,col.names = T,row.names = F,sep="\t")

---

在样本距离聚类时我们发现一些正常组织和肿瘤组织聚集在了一起,根据样本信息,我们发现我们得到的43个正常组织是肿瘤病人的癌旁正常组织,这可能是正常和肿瘤组织聚在一起的原因.我们还发现样本之间的差异过大,生物变异系数过高,这要求我们必须对这些样本进行筛选才能进一步分析.不过我们在阅读文献和查询数据库时发现几乎没有人对样本进行严格过滤,直接拿来分析得出差异表达miRNA结果,我们在后面将过滤后与没过滤的样本差异表达分析的结果进行比对.看是否存在显著差异.

更多原创精彩内容敬请关注生信杂谈：