美文网首页
Reanalysis of mouse ENCODE mRNA-

Reanalysis of mouse ENCODE mRNA-

作者: 宇天_ngs | 来源:发表于2017-10-01 20:20 被阅读68次

前言

ENCODE项目团队发表在PNAS上的文章表示：通过数据分析基因的表达水平在同一个物种不同组织间的相似性比不同物种的相同组织还要高。这与之前的观念完全相悖，不同物种里主要的生物过程、通路等应该是高度保守的，其所具有的相同或类似功能的组织中基因的表达模式应该也是高度相似的。作者分析发现这些数据存在一定的设计缺陷，数据里包含有明显的批效应；通过重分析去除批效应，数据反映出了正确的结论。

A reanalysis of mouse ENCODE comparative gene expression data

注释数据准备

mouse genome: ENSEMBL GRCm38.68.dna_sm.toplevel.fa.gz

mouse GTF: ENSEMBL GRCm38.68.gtf

human genome: iGenomes GRCh37
ENSEMBL官方的human基因组包含了haplotypic region(不属于primary assembly)，比对不考虑这些区域，因此作者选择在iGenomes上下载

human GTF: GENCODE gencode.v14.annotation.gtf

Ortholog annotation(直系同源注释)：http://compbio.mit.edu/modencode/orthologs/modencode.common.orth.txt.gz

批信息提取

通过分析read的标识行，read来自于CASAVA v1.8之后，提取 : 分割的前四列作为批信息，即：machine identifer, run number, flowcell identifer, and flowcell lane number。

共有5批数据，只有1批同时包含了人和鼠的数据。另外4批可分为两对，各具有几乎相同的组织来源。![1506694797631]


1506694797631.png

预处理

通过FastQC发现存在GC偏倚，把overrepresented sequences进行了BLAT，大多都完美匹配到线粒体基因组或者rRNA repeats。也就表明没有污染、接头残留的问题存在。

比对

工具：Tophat v2.0.11

参数：--mate-inner-dist 200 #(PE 100, insert size 350-450nt)
​ --bowtie-n

​ --library-type fr-frststrand #TruSeq Stranded mRNA LT Sample Prep Kit

​ -g 1 #排除多重比对的reads

​ -m 1 #剪接比对时锚定区仅允许1个错配

FPKM

工具：Cufflinks v2.2.1

参数：--library-type fr-frststrand
​ -G

compute fragments per kilo base of transcript per million (FPKM) values and aggregate them per gene. (即先对转录本计算FPKM，再按gene汇总)

同一个名字可能对应多个ID，其中有些是被废除了，针对其余的ID，把表达值合并入首次出现的ID。

FPKM表达矩阵采用log2(FPKM+1)的形式。

Count

工具： featureCounts v1.4.4

参数： -p #计算fragment而不是read
​ -C #嵌合片段不汇总
​ -s 2 #双端reads都是反义链

表达分析

样品热图聚类

suppressMessages(library(pheatmap))
suppressMessages(library(ggplot2))
setwd("C:/Users/Administrator/Desktop/f1000/")
inputDir = "./input"
outputDir = "./output"

print("Performing analysis of FPKM matrix")
#load input
datasets = as.data.frame(scan(paste(inputDir,"Stanford_datasets.txt",sep="/"),list(setname="",seqBatch="",species="",tissue=""),sep="\t"))
fpkmMat <- as.matrix(read.table(paste(inputDir,"Stanford_datasets_fpkmMat.txt",sep="/"),header=FALSE,sep='\t'))

# perform log transformation
logTransformed_fpkmMat = log2(fpkmMat+1)
colnames(logTransformed_fpkmMat) <- datasets$setname

#plot correlation heatmaps
outFname = paste(outputDir,"logTransformed_fpkmMat_pearsonCorrelation_euclideanDistance_completeLinkage.png",sep="/")
png(outFname,width=480,height=480,units="px",pointsize=12)
pheatmap(cor(logTransformed_fpkmMat)) 
#pearson correlation; complete linkage
invisible(dev.off())
# clustering_method:
# "average" average linkage
# "single"  single linkage
# "manhattan" manhattan distance, complete linkage
# "canberra"  canberra distance, complete linkage

pheatmap(cor(logTransformed_fpkmMat,method="spearman"))
# spearman correlation; complete linkage

pheatmap(logTransformed_fpkmMat, clustering_distance_cols="correlation", show_rownames = FALSE) 
# clustering of log-transformed FPKM, Pearson correlation distance, complete linkage
1506772984830.png

可以看出确实都是按物种聚类在了一起。

PCA分析

#pca
transposeLogTransformed_fpkmMat = t(logTransformed_fpkmMat)
# remove genes of zero variance
cat('Performing PCA on log-transformed FPKM matrix',sep='\n')
pca_proc <- prcomp(transposeLogTransformed_fpkmMat[,apply(transposeLogTransformed_fpkmMat, 2, var, na.rm=TRUE) != 0],scale=TRUE,center=TRUE,retx=TRUE)
PoV <- as.numeric(summary(pca_proc)[[6]][2,])
#PoV <- pca_proc$sdev^2/sum(pca_proc$sdev^2)

plotData = datasets[,c("setname","species","tissue")]
plotData$PC1 <- pca_proc$x[,1]
plotData$PC2 <- pca_proc$x[,2]
plotData$PC3 <- pca_proc$x[,3]
plotData$PC4 <- pca_proc$x[,4]
plotData$PC5 <- pca_proc$x[,5]
plotData$PC6 <- pca_proc$x[,6]

# 2D plots of pairs of principal components
i=1; j=2
outFname = paste(outputDir,"/logTransformedFPKM_pca_PC",i,"_vs_PC",j,".png",sep="")
xlabel = paste("PC",i," (",100*round(PoV[i],2),"% variability)",sep="")
ylabel = paste("PC",j," (",100*round(PoV[j],2),"% variability)",sep="")
p=qplot(PC1,PC2,data=plotData,color=species,shape=tissue,xlab=xlabel,ylab=ylabel)+scale_shape_manual(values=c("adipose"=0,"adrenal"=1,"brain"=2,"sigmoid"=3,"heart"=4,"kidney"=5,"liver"=6,"lung"=7,"ovary"=8,"pancreas"=9,"small_bowel"=10,"spleen"=11,"testis"=12))
ggsave(outFname,p,unit="in",width=7,height=5,dpi=600)
logTransformedFPKM_pca_PC1_vs_PC2.png

PCA同样看出人和小鼠物种间的差异很大。

# testing for significance of spearman/pearson correlations between the matched tissues PC values of human and mouse
humanIdx = c(1:nrow(datasets))[datasets$species=="human"]
mouseIdx = c(1:nrow(datasets))[datasets$species=="mouse"]
print("Testing for significance of Spearman correlation between human and mouse PC1 values, matched by tissue:")
cor.test(plotData$PC1[humanIdx],plotData$PC1[mouseIdx],method="spearman")

print("Testing for significance of Pearson correlation between human and mouse PC1 values, matched by tissue:")
cor.test(plotData$PC1[humanIdx],plotData$PC1[mouseIdx],method="pearson")

Raw Counts重分析

suppressMessages(library(pheatmap))
suppressMessages(library(ggplot2))
suppressMessages(library(EDASeq))
suppressMessages(library(edgeR))
suppressMessages(library(sva))

rm(list = ls())
inputDir = "./input"
outputDir = "./output"

#1. load data，数据载入
datasets = as.data.frame(scan(paste(inputDir,"Stanford_datasets.txt",sep="/"),list(setname="",seqBatch="",species="",tissue=""),sep="\t"))
rawCounts <- as.matrix(read.table(paste(inputDir,'Stanford_datasets_rawCountsMat.txt',sep="/"),header=FALSE,sep='\t'))
colnames(rawCounts) <- datasets$setname
geneDetails <- as.data.frame(scan(paste(inputDir,"ortholog_GC_table.txt",sep="/"),skip=1,list(mouse_name="",mouse_GC = 0.0,human_name = "",human_GC=0.0)))
rownames(rawCounts) <- geneDetails$human_name 
rownames(datasets) <- datasets$setname
rownames(geneDetails) <- geneDetails$human_name 

# 先过滤掉表达值位于下30%的基因，再去除线粒体基因
#2. filter out lower 30% + mitochondrial genes
# rowSums = apply(rawCounts,1,function(x) sum(x))
rowSums = rowSums(rawCounts)
low30 <- quantile(rowSums,probs = 0.3) # 2947.9 
filter <- apply(rawCounts,1,function(x) sum(x)>low30 )
mt <- grep("mt-",geneDetails$mouse_name)
filteredNames <- setdiff(rownames(rawCounts[filter,]),rownames(rawCounts[mt,])) 
filteredRawCounts <- rawCounts[filteredNames,]

对过滤后的数据进行GC含量标准化及TMM标准化：

# 调用EDASeq的函数，根据GC含量进行lane间的标准化
#3. Use EDASeq to normalize data within lanes, accounting for GC content
GCnormCounts <- filteredRawCounts
GCnormCounts[,1:13] <- withinLaneNormalization(filteredRawCounts[,1:13],geneDetails[filteredNames,"human_GC"],which="loess",round=TRUE)
GCnormCounts[,14:26] <- withinLaneNormalization(filteredRawCounts[,14:26],geneDetails[filteredNames,"mouse_GC"],which="loess",round=TRUE)

# 调用edgeR的TMM方法计算测序深度归一化因子
#4. depth normalize,using TMM scaling factors - divide by sum, then multiply by mean of sums
# origColSums <- apply(rawCounts,2,function(x) sum(x))
origColSums <- colSums(rawCounts)
normFactors <- calcNormFactors(GCnormCounts,method='TMM')
normalizedColSums <- origColSums
i <- 1
while (i<length(colSums)){
  normalizedColSums[i] <- origColSums[i]* normFactors[i]
  i <- i+1
}
# ??? 最后一个样本不标准化 ???
# normalizedColSums <- origColSums*normFactors

meanDepth <- mean(normalizedColSums)
filteredDepthNormCounts <- GCnormCounts
i <- 1
while (i<ncol(filteredDepthNormCounts)){
  filteredDepthNormCounts[,i] <- (GCnormCounts[,i]/normalizedColSums[i])*meanDepth
  i <- i+1
}

filteredDepthNormCounts <- t(apply(GCnormCounts,1,function(x) (x/normalizedColSums)*meanDepth))
                                   
#5. log transformation
logTransformedDepthNormCounts <- log2(filteredDepthNormCounts+1)

数据用sva移除批效应：

#6. use combat on log2 values to remove batch effects 
meta <- data.frame(seqBatch = datasets$seqBatch,tissue=datasets$tissue,species=datasets$species)
design <- model.matrix(~1,data=meta)
combat <- ComBat(dat = logTransformedDepthNormCounts,batch=meta$seqBatch,mod=design,par.prior=TRUE)

批效应移除后再通过PCA和热图查看效果：

#7.plot  heatmaps
outFname = paste(outputDir,"batchCorrected_logTransformedProcessedCounts_pearsonCorrelation_euclideanDistance_completeLinkage.png",sep="/")
png(outFname,width=480,height=480,units="px",pointsize=12)
pheatmap(cor(combat)) #pearson correlation; complete linkage; euclidean distance
invisible(dev.off())

#8. PCA
transposed_combat <- t(combat)
pca_proc <- prcomp(transposed_combat[,apply(transposed_combat, 2, var, na.rm=TRUE) != 0],scale=TRUE,center=TRUE,retx=TRUE)
PoV <- as.numeric(summary(pca_proc)[[6]][2,])

plotData = datasets[,c("setname","species","tissue")]
plotData$PC1 <- pca_proc$x[,1]
plotData$PC2 <- pca_proc$x[,2]
plotData$PC3 <- pca_proc$x[,3]
plotData$PC4 <- pca_proc$x[,4]
plotData$PC5 <- pca_proc$x[,5]
plotData$PC6 <- pca_proc$x[,6]

# 2D plots of pairs of principal components
i=1; j=2
outFname = paste(outputDir,"/batchCorrected_logTransformedProcessedCounts_pca_PC",i,"_vs_PC",j,".png",sep="")
xlabel = paste("PC",i," (",100*round(PoV[i],2),"% variability)",sep="")
ylabel = paste("PC",j," (",100*round(PoV[j],2),"% variability)",sep="")
p=qplot(PC1,PC2,data=plotData,color=species,shape=tissue,xlab=xlabel,ylab=ylabel)+scale_shape_manual(values=c("adipose"=0,"adrenal"=1,"brain"=2,"sigmoid"=3,"heart"=4,"kidney"=5,"liver"=6,"lung"=7,"ovary"=8,"pancreas"=9,"small_bowel"=10,"spleen"=11,"testis"=12))
ggsave(outFname,p,unit="in",width=7,height=5,dpi=150)
1506847646436.png

移除批效应后可以发现人和鼠相同组织间的相似度更高了！

1506847753917.png

PCA图中样品间的物种差异度更小(混合程度更高)！

PS：人的样品因收集过程繁琐，其RIN值可能更低，这样就需要根据GC含量进行标准化，再进行其他标准化；另外在本实验设计中，批效应和物种带来的差异完全混杂，无法区分，因为这种去除批效应的方式还是有争议的。

相关文章

网友评论

      本文标题：Reanalysis of mouse ENCODE mRNA-

      本文链接：https://www.haomeiwen.com/subject/mqvaextx.html

      延伸阅读

      深度阅读

      • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

      栏目导航

      热点阅读

      关于我们|服务条款|联系我们|Reanalysis of mouse ENCODE mRNA-|投稿指南|网站地图|RSS订阅|排版工具|手机版

      提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

      Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

      备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

      本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

      所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！