前言
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含量进行标准化，再进行其他标准化；另外在本实验设计中，批效应和物种带来的差异完全混杂，无法区分，因为这种去除批效应的方式还是有争议的。
