GTEX表達(dá)矩陣預(yù)處理

#設(shè)置工作路徑
setwd("D:/gtex")
#初始化環(huán)境：
rm(list=ls())
#載入包：
library(data.table)
#獲得要提取的樣本ID：
sample_id <- read.table("sample.txt")
#再添加兩個(gè)列名：
sample_id <-append(c("Name","Description"),sample_id)
#1.3G的表達(dá)譜文件路徑：
GTExFile_path<-"./GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_reads.gct"
#根據(jù)要提取的樣本ID讀入文件,skip前兩行，select要讀取的列，將第一行用作行名，分隔符Tab
GTEx_reads<-fread(GTExFile_path,sep="\t",header =T,skip = 2,select = sample_id)
#儲(chǔ)存ensembl基因名與gene symbol的對(duì)應(yīng)關(guān)系
#write.table(GTEx_reads[,1:2],file = "vs.tsv")
#刪除前兩列數(shù)據(jù)
GTEx_reads_ <- GTEx_reads[,-(1:2)]
#將字符改為數(shù)字
GTEx_reads_1 <- as.data.frame(lapply(GTEx_reads_,as.numeric))
#用ENSEMBL基因名作行名
vs <- read.table("vs.tsv")
rownames(GTEx_reads_1) <- vs[,1]
#更改列名，將字符改為數(shù)字不小心也更改了列名
colnames(GTEx_reads_1) <- colnames(GTEx_reads_)
#將數(shù)據(jù)保存為Rdata
save(GTEx_reads_1,file = "step_1.Rdata")

R語言學(xué)習(xí)：
1.append添加新的列。
2.R語言處理大規(guī)模數(shù)據(jù)速度不算快，通過安裝其他包比如data.table可以提升讀取處理速度。fread比read.table快4—5倍數(shù)。

原始數(shù)據(jù)

GTEX_reads

利用edgeR包對(duì)數(shù)據(jù)進(jìn)行過濾以及歸一化

#設(shè)置工作路徑
setwd("D:/gtex")
#檢查內(nèi)存分配限制
memory.limit()
#設(shè)置為1G內(nèi)存
memory.limit(1000000)
#初始化環(huán)境：
rm(list=ls())
#安裝R包
if(F)
{
  if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
  BiocManager::install("edgeR")
}
library(edgeR)
#導(dǎo)入數(shù)據(jù)，創(chuàng)建DGElist
load(file = "step_1.Rdata")
x <- GTEx_reads_1
#x <- x[,1:560]
tmp <- read.csv(file = "sample_id&tissue.csv")
tmp_1 <- tmp[tmp[,1] %in% colnames(x),]
group = factor(tmp_1[,2])
y = DGEList(counts = x,group = group)
dim(y)
#手動(dòng)過濾
#table(rowSums(y$counts==0)==1)
cut.off.cpm = 10*10**6/median(colSums(x))
 keep_cpm <- rowSums(cpm(y)>cut.off.cpm) >= 36
y_1 <- y[keep_cpm,, keep.lib.sizes=FALSE]
dim(y_1)
#自動(dòng)過濾
if(F){
  keep.exprs <- filterByExpr(y, group=group)
  y_2 <- y[keep.exprs,, keep.lib.sizes=FALSE]
  dim(y_2)
}

y_3 <- calcNormFactors(y_1, method = "TMM")
#y_1$samples$norm.factors

#手動(dòng)得到有效庫
library("dplyr")
y_3[["samples"]] %>% mutate(lib.size=lib.size*norm.factors) -> y_3[["samples"]]

#計(jì)算過濾以及歸一化之后的cpm值
cpm_3 <- cpm(y_3)
if(F)
{
  #voom從x中自動(dòng)提取文庫大小和歸一化因子，以此將原始計(jì)數(shù)轉(zhuǎn)換為log-CPM值。
  design = model.matrix(~group)
  v <- voom(y_3, design)
  #比較
  lcpm <- cpm(y_3, log=TRUE)
  #相似但是不完全相等，因?yàn)関oom使用了更小的預(yù)先計(jì)數(shù)
}

#更改列名為組織
table <- as.data.frame(table(y_3[["samples"]][,1]))
b <- character()
for(i in 1:nrow(table) ){
  a <- rep(levels(y_3[["samples"]][,1])[i] , each = table[,"Freq"][i])
  b <- c(b,a)
}
colnames(cpm_3) <- b

#求每個(gè)基因在每個(gè)組織中表達(dá)的中位值
library(reshape2)
melt <- melt(cpm_3)
#aggregate(melt,by=list(melt$Var1,melt$Var2),FUN=median)
by <- group_by(melt,Var1,Var2) %>% summarise(median=median(value))
#數(shù)據(jù)長寬轉(zhuǎn)換
library(tidyfst)
by %>% wider_dt("Var1",name = "Var2",value = "median") -> melt_1
save(melt_1,file = "step_2.Rdata")

樣本信息：

> table(tmp_1[,2])

Adrenal Gland         Brain        Breast         Heart        Kidney         Liver 
          103           291           166           207            36           100 
         Lung        Muscle      Pancreas          Skin       Stomach        Testis 
          208           367           125           597           112           175 
       Uterus 
           48 

一共得到了13個(gè)組織，2535個(gè)樣本。
低表達(dá)基因過濾：

> dim(y_1)
[1] 36692  2535

得到了36692個(gè)基因。

幾點(diǎn)說明：
1.手動(dòng)過濾的閾值選擇：
此函數(shù)默認(rèn)選取最小的組內(nèi)的樣本數(shù)量為最小樣本數(shù)，保留至少在這個(gè)數(shù)量的樣本中有10個(gè)或更多計(jì)數(shù)的基因。實(shí)際進(jìn)行過濾時(shí)，使用的是CPM值而不是表達(dá)計(jì)數(shù)，以避免對(duì)總序列數(shù)大的樣本的偏向性。在這個(gè)數(shù)據(jù)集中，總序列數(shù)的中位數(shù)是5.1千萬，且10/51約等于0.2，所以filterByExpr函數(shù)保留在至少三個(gè)樣本中CPM值大于等于0.2的基因。在我們的此次實(shí)驗(yàn)中，一個(gè)具有生物學(xué)意義的基因需要在至少三個(gè)樣本中表達(dá)，因?yàn)槿N細(xì)胞類型組內(nèi)各有三個(gè)重復(fù)。過濾的閾值取決于測序深度和實(shí)驗(yàn)設(shè)計(jì)。如果樣本總表達(dá)計(jì)數(shù)量增大，那么可以選擇更低的CPM閾值，因?yàn)楦蟮目偙磉_(dá)計(jì)數(shù)量提供了更好的分辨率來探究更多表達(dá)水平更低的基因。
對(duì)于在大多數(shù)樣本中表達(dá)數(shù)量都很少的基因（往往是沒有生物學(xué)意義的基因），需要進(jìn)行過濾，這一步可以根據(jù)自己定義的標(biāo)準(zhǔn)過濾，edgeR推薦使用該包的CPM（ count-per-million ）值進(jìn)行過濾。
在這里我也使用了cpm值進(jìn)行手動(dòng)的過濾。

2.歸一化方法：
我們需要使用 edgeR 中的 calcNormFactors 函數(shù)，用 TMM(Robinson and Oshlack 2010) 方法進(jìn)行歸一化。此處計(jì)算得到的歸一化系數(shù)被用作文庫大小的縮放系數(shù)。當(dāng)我們使用DGEList對(duì)象時(shí)，這些歸一化系數(shù)被自動(dòng)存儲(chǔ)在xnorm.factors。對(duì)此數(shù)據(jù)集而言，TMM歸一化的作用比較溫和，這體現(xiàn)在所有的縮放因子都相對(duì)接近1。這種歸一化的方法可以消除RNA組成所帶來的誤差。
(TMM方法文獻(xiàn)：https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-3-r25）

標(biāo)準(zhǔn)化因子的生物學(xué)意義
其實(shí)這個(gè)標(biāo)準(zhǔn)化因子算法就是選出一個(gè)有代表性的gene X（其實(shí)是每個(gè)樣本一個(gè)代表性gene X），而這個(gè)gene X的reads for gene X/average for gene X比值就是標(biāo)準(zhǔn)化因子。
只不過選取gene X的時(shí)候，通過對(duì)數(shù)變換和中位數(shù)的方法，更多的參考了中等表達(dá)基因和管家基因的數(shù)據(jù)趨勢，而剔除了特異性表達(dá)基因和高差異表達(dá)基因的影響。
相比較RPKM,FPKM,TPM標(biāo)準(zhǔn)化方法是除以總Read數(shù)，DESeq2標(biāo)準(zhǔn)化方法是除以一個(gè)有代表性基因的Read數(shù)，只不過這個(gè)Read數(shù)進(jìn)行了變換（它除以了幾何平均Read數(shù)， reads for gene X/average for gene X）。因?yàn)楦芴幚泶嬖谔禺愋员磉_(dá)基因和高差異表達(dá)基因的數(shù)據(jù)。

通過得到calcNormFactors，我們可以得到有效庫大小，從而得到消除了測序深度以及在某個(gè)樣本中特異性表達(dá)基因（RNA組成）所帶來的影響。

3.最終使用CPM值做接下來的分析。選取樣本的中位值做分析。（中位值參考文獻(xiàn)。）

4.R語言學(xué)習(xí)：tidyfst、reshape2包進(jìn)行數(shù)據(jù)的長寬轉(zhuǎn)換，dplyr包中的mutate、group_by函數(shù)進(jìn)行列的計(jì)算、分組計(jì)算等。

整理為規(guī)范的SPM上游分析文件

#設(shè)置工作路徑
setwd("D:/gtex")
#檢查內(nèi)存分配限制
memory.limit()
stringsAsFactors=FALSE
#設(shè)置為1G內(nèi)存
memory.limit(1000000)
#初始化環(huán)境：
rm(list=ls())
#導(dǎo)入數(shù)據(jù)
vs <- read.table(file = "vs.tsv")
load(file = "step_2.Rdata")
#將基因名導(dǎo)入數(shù)據(jù)框中并用作行名，多個(gè)ensembl_id對(duì)應(yīng)一個(gè)基因名就取其中一個(gè)
Var1 = as.character(melt_1$Var1)
melt_1 <- melt_1[,-1]
melt_1 <- data.frame(cbind(Var1,melt_1))
#tmp <- vs[vs[,1] %in% melt_1[,1],c(1:2)]
colnames(melt_1)[1] <- "Name"
#tmp_1 <- merge(tmp,melt_1,by = "Name")
tmp_1 <- merge(vs,melt_1,by = "Name")
length(unique(tmp_1[,2]))
library(dplyr)
tmp_2 <- distinct(tmp_1,Description,.keep_all = TRUE)
rownames(tmp_2) <-  tmp_2[,2]
exprset <- tmp_2[,-c(1:2)] 
#刪除中間文件
rm(tmp,tmp_1,tmp_2)
colnames(exprset) <- c(
  "Adrenal_Gland",
  "Brain",
  "Breast",
  "Left_Ventricle",
  "Kidney",
  "Liver",
  "Lung",
  "Skeletal_Muscle",
  "Pancreas",
  "Skin",
  "Stomach",
  "Testis",
  "Uterus"
)

write.table(exprset,file = "exprset.txt",sep = '\t')

> dim(exprset)
[1] 36598    13

說明：如果多個(gè)ensembl name對(duì)應(yīng)一個(gè)gene symbol，我們就只取其中的一個(gè)ensembl name進(jìn)行后續(xù)的分析。這一步的基因數(shù)目從36692到36598，誤差不大。
warning：這一步涉及到手動(dòng)改列名，發(fā)現(xiàn)了一個(gè)錯(cuò)誤，以后一定要多加小心，能不手動(dòng)盡量不要手動(dòng)。