之前小云和小伙伴分享了如何使用TCGA數(shù)據(jù)庫(kù)，同時(shí)以胰腺癌數(shù)據(jù)為實(shí)例進(jìn)行了數(shù)據(jù)下載，那么接下來(lái)就要對(duì)下載的數(shù)據(jù)根據(jù)需要進(jìn)行整理

TCGA數(shù)據(jù)整理

# 數(shù)據(jù)讀取 -------------------------------------------------------------------

首先讀取下載數(shù)據(jù)的記錄，這樣免得大家二次下載，省去等待時(shí)間

load("TCGA_PAAD.Rdata")

?

# 加載所需要的包 -------------------------------------------------------------

library(dplyr) #數(shù)據(jù)整理

#只取protein coding，編碼蛋白的mRNA進(jìn)行分析

table(expr_count$gene_type)

mRNA <- expr_count[expr_count$gene_type=="protein_coding",]

?

# count數(shù)據(jù)清洗 ---------------------------------------------------------------

mRNA <- mRNA[,-1]

dim(mRNA)#查看一共多少行

mRNA$mean <- rowMeans(mRNA[,2:184])#計(jì)算每一行表達(dá)量的平均值

mRNA <- arrange(mRNA,desc(mean))#按mean值大小降序排列

mRNA <- mRNA %>% distinct(gene_name, .keep_all = T)#刪除重復(fù)Gene的行

mRNA <- select(mRNA,-c(mean))#刪除之前生成的mean一行

?

#將腫瘤和非腫瘤樣本區(qū)分開排序

rownames(mRNA) <- mRNA$gene_name

mRNA <- mRNA[,-1]

table(substr(colnames(mRNA),14,16))

那么-01A和-11A，這里的字母A又是什么含義呢？

從這個(gè)示意圖上可以看到，有時(shí)候我們可以從一個(gè)病例身上取多個(gè)樣本，不論是腫瘤樣本，還是癌旁正常對(duì)照，然后存放在不同的管子里面。這里的A，B，C就表示樣本的順序。官方文檔的解釋如下。

我們可以看到后綴中數(shù)字與樣本類型的對(duì)應(yīng)關(guān)系，如下圖

解釋完樓，來(lái)我們繼續(xù)~

Tumor <- grep('01A',colnames(mRNA))

Tumor ?#腫瘤樣本所處位置

Tumor_mRNA <- mRNA[,Tumor]

?

Normal <- grep('11A',colnames(mRNA))

Normal #正常樣本所處位置

Normal_mRNA <- mRNA[,Normal]

?

mRNA_count <- cbind(Normal_mRNA,Tumor_mRNA)

write.csv(mRNA_count,file = "mRNA_count.csv")

?

# TPM數(shù)據(jù)清洗 -----------------------------------------------------------------

mRNA_TPM <- expr_TPM[expr_TPM$gene_type=="protein_coding",]

mRNA_TPM <- mRNA_TPM[,-1]

dim(mRNA_TPM)#查看一共多少行

mRNA_TPM <- mRNA_TPM %>%

??mutate(mean=rowMeans(.[,-1])) %>%

??arrange(desc(mean)) %>%

??distinct(gene_name,.keep_all = T) %>%

??select(-mean)

#將腫瘤樣本和非腫瘤樣本分開排序

rownames(mRNA_TPM) <- mRNA_TPM$gene_name

mRNA_TPM <- mRNA_TPM[,-1]

table(substr(colnames(mRNA_TPM),14,16))

?

Tumor <- grep('01A',colnames(mRNA_TPM))

Tumor #腫瘤樣本所處位置

Tumor_mRNA_TPM <- mRNA_TPM[,Tumor]

?

Normal <- grep('11A',colnames(mRNA_TPM))

Normal #正常樣本所處位置

Normal_mRNA_TPM <- mRNA_TPM[,Normal]

?

mRNA_TPM_new <- cbind(Normal_mRNA_TPM,Tumor_mRNA_TPM)

?

# 保存文件用于下次的分析 -------------------------------------------------------------

save(mRNA_count,mRNA_TPM_new,file = "count_and_TPM.Rdata")

?

DESeq2差異分析

數(shù)據(jù)整理完成后，一般我們會(huì)獲得一個(gè)原始?read count表達(dá)矩陣，其中行是基因，列是樣品。常用的差異分析工具包括limma，edgeR，DESeq2。小圖就以R包DESeq2的差異表達(dá)基因分析方法為例做個(gè)簡(jiǎn)單演示，因?yàn)樾D今天下載的是count數(shù)據(jù)。

簡(jiǎn)單地說(shuō)，DESeq2將對(duì)原始reads count進(jìn)行建模，使用標(biāo)準(zhǔn)化因子（scale factor）來(lái)解釋庫(kù)深度的差異。然后DESeq2估計(jì)基因的離散度，并縮小這些估計(jì)值以生成更準(zhǔn)確的離散度估計(jì)，從而對(duì)reads count進(jìn)行建模。最后，DESeq2擬合負(fù)二項(xiàng)分布的模型，并使用Wald檢驗(yàn)或似然比檢驗(yàn)進(jìn)行假設(shè)檢驗(yàn)。

來(lái)，馬上開始實(shí)操~

setwd("D:/ DESEq2差異表達(dá)")

?

# 清空變量 --------------------------------------------------------------------

rm(list=ls(all=TRUE))

options(stringsAsFactors = F)

library(DESeq2)

?

expr <- expr[rowMeans(expr) > 1,]

rownames(expr)=mRNA_count[,1]

# 創(chuàng)建分組信息 ------------------------------------------------------------------

table(substr(colnames(expr),14,16))#查看腫瘤樣本和正常樣本數(shù)量

Tumor <- grep('01A',colnames(expr))#查看腫瘤樣本所處位置

Tumor

group <- factor(c(rep("Normal",times=4),rep("Tumor",times=178)))#創(chuàng)建分組因子變量

?

Data <- data.frame(row.names = colnames(expr), #創(chuàng)建分組數(shù)據(jù)框

???????????????????group = group)

View(Data)

#基因表達(dá)矩陣(expr)、樣本分組數(shù)據(jù)框(Data)都已經(jīng)準(zhǔn)備完畢

?

?

# 開始進(jìn)行差異表達(dá)分析 --------------------------------------------------------------

#第一步：構(gòu)建DEseq2對(duì)象(dds)

dds <- DESeqDataSetFromMatrix(countData = expr,

??????????????????????????????colData = Data,

??????????????????????????????design = ~ group)

#第二步：開始差異分析

dds2 <- DESeq(dds)

res <- results(dds2, contrast=c("group", "Tumor", "Normal"))#腫瘤在前，對(duì)照在后

##或者res= results(dds)

res <- res[order(res$pvalue),]

summary(res)

my_result <- as.data.frame(res)#轉(zhuǎn)成容易查看的數(shù)據(jù)框

my_result <- na.omit(my_result)#刪除倍數(shù)為0的值

#第三步：保存差異分析的結(jié)果

library(dplyr)

my_result$Gene_symbol<-rownames(my_result)

my_result <- my_result %>% dplyr::select('Gene_symbol',colnames(my_result)[1:dim(my_result)[2]-1],everything())

rownames(my_result) <- NULL

write.csv(my_result,file="my_result_DESeq2.csv")

?

?

# DEG的篩選 ------------------------------------------------------------------

my_result$regulate <- ifelse(my_result$pvalue > 0.05, "unchanged",

???????????????????????ifelse(my_result$log2FoldChange > 0, "up-regulated",

??????????????????????????????ifelse(my_result$log2FoldChange < 0, "down-regulated", "unchanged")))

table(my_result$regulate)

#可以把上調(diào)基因和下調(diào)基因取出放在一塊

DEG_deseq2 <-subset(my_result, pvalue < 0.05 & abs(log2FoldChange) > 0)

upgene <- DEG_deseq2[DEG_deseq2$regulate=='up-regulated',]

downgene <- DEG_deseq2[DEG_deseq2$regulate=='down-regulated',]

write.csv(DEG_deseq2,file= "DEG_deseq2_0.05.csv")

好了今天的進(jìn)階版TCGA數(shù)據(jù)整理分享就到這啦，小伙伴們有遇到問(wèn)題歡迎隨時(shí)來(lái)找小云分享討論呀！