R根据基因的表达量筛选基因

data <- read.delim('mouse_gc_log2.xls', row.names = 1, sep = '\t', check.names = FALSE)#行为基因ID列为样本的表达矩阵
data[order(rowSums(data),decreasing=T)[1:3000],]  #筛选表达量高的前3000个基因
data[rowSums(data)>1,]  #过滤掉表达量低的基因

data1=rowMeans(data)>1  #按平均数筛选
data2=rowSums(data>0)>6 #表达量不为0的样品个数筛选
data=data[data1 & data2,] #联合一下
data$cv <- apply(data, 1, function(x){
  sd(x)/mean(x)*100
})
data_df <- data[order(data$cv, decreasing = T)[1:3000], 1:15]#筛选变异系数最大的3000个基因

Gene-rank score计算

data <- read.delim('mouse_gc.xls', row.names = 1, sep = '\t', check.names = FALSE)
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b <- apply(data, 2, function(y) rank(y) / length(y))
write.table(b, file="gene_rank_score.xls", sep="\t", quote=F, row.names=T, col.names=T)
#gene-rank没有改变样本间的相关性

MetabolAnalyze实现Pareto scaling

BiocManager::install("MetabolAnalyze")
library("MetabolAnalyze")
data <- read.delim('mouse_gc_TMM_cross.xls', row.names = 1, sep = '\t', check.names = FALSE)
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b = scaling(data,type = "pareto")
boxplot(b)
write.table(b, file="Pareto_scaled.xls", sep="\t", quote=F, row.names=T, col.names=T)

edgeR之TMM标准化

library(edgeR)
library(ggplot2)
data <- read.delim('../mouse_TPM_notCross.matrix.xls', row.names = 1, sep = '\t', check.names = FALSE)
group <- factor(rep(c('A', 'B', 'C'), each = 3))
y <- DGEList(counts = data, group = group)
y <- calcNormFactors(y)
logcpm <- cpm(y, prior.count=3, log=TRUE)
write.table(logcpm, file="mouse_TMM.xls", sep="\t", quote=F, row.names=T, col.names=T)
#dgList <- estimateCommonDisp(y)
#dgList <- estimateTagwiseDisp(dgList)
#norm_counts.table <- t(t(dgList$pseudo.counts)*(dgList$samples$norm.factors))
#write.table(norm_counts.table, file="mouse_gc_pareto_TMM.xls", sep="\t", quote=F)
#pheatmap(log(data+1),cluster_rows=T,cluster_cols=T,scale="none",border_color="white",color=colorRampPalette(rev(c("red","white","blue")))(102))

limma去除批次效应

什么是批次效应(batch effect)?
不同平台的数据,同一平台的不同时期的数据,同一个样品不同试剂的数据,以及同一个样品不同时间的数据等等都会产生一种batch effect 。这种影响如果广泛存在应该被足够重视,否则会导致整个实验和最终的结论失败。比对实验组和对照组,不同的处理是患病和不患病(测序时,先测得疾病,然后测得正常),然后你通过分析,得到很多差异表达的基因。现在问题来了,这个差异表达的结果是和你要研究的因素有关,还是时间有关,这个问题里时间就会成为干扰实验结果的因素,这个效应就是batch effect。

library(limma)
data <- read.table("mouse_gc.xls",sep="\t",header = TRUE)
batch1 <- c(rep('fish',9),rep('mouse',9))
batch1 <- as.factor(batch1)
design <- model.matrix(~0 + batch1)
#data=normalizeBetweenArrays(data)#如果组内的中位数不在同一条水平线上,现在该参数校正,再使用下面的批次校正 new_data <- removeBatchEffect(data[,2:19], batch = batch1) boxplot(data[,2:19]) boxplot(new_data) write.table(new_data, file="removeBatch.xls", sep="\t", quote=F, row.names=T, col.names=T)

limma结果出现负值请参考:https://cloud.tencent.com/developer/article/1680305
参考链接:https://www.omicsclass.com/article/1113

preprocessCore分位数标准化

library("preprocessCore")
library("ggplot2")
data <- read.table("mouse_gc.xls",sep="\t",header = TRUE)
rownames(data) <- data$GeneID
#rownames(data) <- data[,1]
#data_mat <- data.matrix(data[,-1]) 
quantile_normalisation <- function(df){
  df_rank <- apply(df,2,rank,ties.method="min")
  df_sorted <- data.frame(apply(df, 2, sort))
  df_mean <- apply(df_sorted, 1, mean)

  index_to_mean <- function(my_index, my_mean){
    return(my_mean[my_index])
  }

  df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean)
  rownames(df_final) <- rownames(df)
  return(df_final)
}
new_data <- quantile_normalisation(data[,2:19])
boxplot(data)
boxplot(new_data)
write.table(new_data, file="111.xls", sep="\t", quote=F, row.names=T, col.names=T)
##############################################################
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b=normalize.quantiles(data)
boxplot(b)
##使用normalize.quantiles函数和上面结果相同

跨物种比较转录组

近几年几篇高分文献对跨物种比较转录组方法的应用及研究领域,感兴趣的可自行查看原文献:

Cell_2019_跨物种单细胞测序揭示小胶质细胞基因表达谱进化特征

Science_2020_人类猪小鼠大脑中的蛋白编码基因图谱

Science_2018_跨人类和猕猴大脑发育的时空转录组差异

Cell_2021_African lungfish genome sheds lighton the vertebrate water-to-landtransition

Cell_2021_人类大脑进化的发育调控机制

Cell_2021_原始辐鳍鱼类基因组的解析发现登陆的遗传基础已经在硬骨鱼类祖先出现

看完上述文献可知和种内比较转录组相比,跨物种(种间)比较转录组要复杂的多,主要体现在背景基因总数、测序深度、管家基因表达谱和可能的其他因素差异等;分析流程主要分为两步:1、直系同源基因的查找。此处大部分文献使用 one-to-one orthologs (及双向最优比对),如果物种跨度较大可选较为宽松的阈值,期望值1e-10,覆盖度50%;如果物种亲缘关系较为接近就选较严的阈值。也有文献先获得直系同源基因列表,一个gene list及为一个基因(matagene),该基因的表达量为gene list所有基因表达量的和(基因拷贝数差异)。2、合并后基因表达矩阵的标准化。标准化方法各异,方法之一:首先在种内计算出FPKM、RPKM、TPM等,种内TMM标准化,然后根据背景基因集合并数据得到新的表达矩阵再使用TMM标准化。也有直接quantile normalization using the R package preprocessCore(具体R实现可参考下篇文章),标准化方法的选择对结果的影响很大,可自行尝试适合自己的标准化方法。另外,不同数据来源需要去除批次效应,不然进行PCA和样本相关性分析时会将不同数据来源的样本聚到一块;数据的适当缩放也是必要的,一些聚类方法会将极值带来的影响放大化,使结果不准确。下图是Science_2020_人类猪小鼠大脑中的蛋白编码基因图谱分析流程。