通過整理TCGA數據，探索某癌症的癌組織和正常組織的差異基因。
實驗設計
實驗目的決定試驗方法和途徑。
試驗目的 :獲取三陰性乳腺癌的正常組織和癌症組織的基因表達差異情況，比較三陰性乳腺癌中的基因表達變化情況。
試驗設計 :通過TCGA獲取乳腺癌的RNA-seq表達數據，篩選出三陰性乳腺癌的樣本，通過比較癌症和正常組織的表達差異。
TCGA資料庫簡介
一句話介紹:TCGA資料庫是一個由國家癌症研究所(National Cancer Institute)和美國人類基因組研究所(National Human Genome Research Institute)共同監督的一個項目。使用對患者樣本的高通量基因組測序和分析技術來試圖提供括基因表達譜，拷貝數變異分析，SNP基因分型，全基因組DNA甲基化分析，微RNA分析等信息。收錄了33種癌症基因組測序數據。TCGA數據處理和整理比Oncoman和GEO困難一些。但是針對腫瘤和癌症所能提供的信息是很完善和可靠的。 TCGA和GEO存在的區別是，GEO存在各種研究領域和研究方向的NGC數據和分析。TCGA是專門針對腫瘤和癌症設立的。TCGA優勢是豐富且規範的臨床數據，以及針對每種癌型的大樣本量。
TCGA數據的獲取
背景介紹： 三陰性乳腺癌是指癌組織免疫組織化學檢查結果為雌激素受體（ER）、孕激素受體（PR）和原癌基因Her-2均為陰性的乳腺癌。這類乳腺癌占所有乳腺癌病理類型的10.0%～20.8%，具有特殊的生物學行為和臨床病理特徵，預後較其他類型差。--from:百度百科:三陰性乳腺癌
實驗設計：
三陰性乳腺癌的篩選標準是根據pheotype來確定的，表達量是通過RNAseq結果確定的，正常組織和癌症組織是通過病例號確定的。 1. TCGA項目的數據可以通過Genomic Data Commons Data Portal獲取，即通過GDC來訪問，訪問地址：https://portal.gdc.cancer.gov/。 2. TCGA資料庫公開免費，所以有許多針對TCGA數據進行整合的網站。還可以通過UCSC Xena進行下載：https://xena.ucsc.edu/public。試驗數據的選擇和目的息息相關，試驗設計：
點擊BRCA（乳腺癌）進入數據選擇界面。 3. 選擇gene expression RNAseq>HTSeq - Counts。注意不要用RPKM等經過了normlization的表達矩陣來分析。要使用Counts來進行差異分析，因為在差異分析時候會自動進行標準化。如果數據經過處理例如log2+1，則可以下載後逆運算轉變回來。 4. 選擇phenotype>Phenotype這裡面有病人的病例信息等，可以通過統計篩選出三陰性乳腺癌的患者。
5. 點擊連接進去會看到詳細的信息如下載地址、樣品數、數據處理方法等。
數據預處理
實驗設計：篩選出三陰性乳腺癌的患者ID，再篩選出同時有癌症組織樣本和癌旁組織樣本，計算初始Count值。
三陰性乳腺癌患者（TNBC）篩選
實驗設計：在R語言中處理數據，選擇breast_carcinoma_estrogen_receptor_status（ER）、PR、HER2受體一欄全為隱形（Negative）的患者。 原始的phenotype文件如下圖，信息量巨大
R語言實現：
讀取文件
phenotype_file <- read.table("TCGA-BRCA.GDC_phenotype.tsv",header = T, sep = , quote = "")
phenotype_colnames <- as.data.frame(colnames(phenotype_file))
#讀取文件並去讀列名（penotype的類型）
table(phenotype_file$breast_carcinoma_estrogen_receptor_status)
#查看雌激素受體(ER)的表達狀態。
table(phenotype_file$breast_carcinoma_progesterone_receptor_status)
#查看孕激素受體（PR）狀態
table(phenotype_file$lab_proc_her2_neu_immunohistochemistry_receptor_status)
#查看原癌基因Her-2狀態
列出三項指標的列表，方便篩選
實驗設計：查找三項指標的列並列出
colnames_num <- grep("receptor_status",colnames(phenotype_file))
#在phenotype_file列中檢索「receptor_status」,grep返回值是列名的列表中包含有「receptor_status」的列序號，因為三陰性乳腺癌的三項指標在這裡都有phenotype_file欄位
#>colnames_num
#[1] 20 25 67 77 87 93
phenotype_colnames <- colnames(phenotype_file)[colname_num]
#利用匹配的返回值取列名到phenotype_colnames中,phenotype_colnames中存儲的是含有「receptor_status」欄位的列名
eph <- phenotype_file[,colnames_num[1:3]]
#把三項指標篩選出來，建立一個新的表格。方便篩選全陰的行（ID號）
函數學習：
?apply
# b為：
# first second
# one 1 4
# two 2 5
# three 3 6
# apply(b,1,sum)
# 上面的指令代表對矩陣b進行行計算，分別對每一行進行求和。函數涉及了三個參數：
# 第一個參數是指要參與計算的矩陣；
# 第二個參數是指按行計算還是按列計算，1——表示按行計算，2——按列計算；
# 第三個參數是指具體的運算參數。
# 上述指令的返回結果為：
# one two three
# 5 7 9
TNBC的篩選
實驗設計：查找指標的狀態並統計；列出三陰性並篩選出來這些數據
tnbc_rownum <- apply(eph,1,function(x)sum(x=="Negative"))
#eph：記錄三項指標類型的矩陣
#將eph按照列執行sum(x=="Negative"):統計各行中，eph列中有Negative的數量
#通過查看陰性的數量，函數返回結果是一個數字向量。內容是0或1、2、3。注意數字的排列順序是對應病人ID的行。
tnbc_rownum
tnbc_sample <- phenotype_file[tnbc_rownum == 3, 1]
#通過讓 tnbc_rownum==3 判別，如果成立，會返回一個logical向量，logical向量是可以直接被矩陣引用的，只會讀取TURE的行或著列
#統計每一行中數值是3的（即全為「Negative」）並取第一列即病人ID
#tnbc_sample：存儲著三陰性乳腺癌的
tnbc_sample
save(tnbc_sample,file = tnbc_sample.Rdata)
#保存Rdata
基因表達矩陣的構建
基因表達矩陣的讀取和讀取後格式修改
實驗設計：轉換data.frame ；行名修改；Count值還原
library(data.table)
#?data.table
#R中的data.table包提供了一個data.frame的高級版本，讓你的程序做數據整型的運算速度大大的增加。
a <- fread("TCGA-BRCA.htseq_counts.tsv.gz",sep = ,header = T)
#讀取表達矩陣文件
a <- as.data.frame(a)
#轉換為數據框，讀取後，轉換前的a的類型是"data.table" "data.frame"
a[1:4,1:4]
rownames(a) <- a[,1]
a <- a[,-1]
#去除第一列，行名修改。
genes <- row.names(a)
genes[1:10]
##接下來還原counts數，網站使用了log2(count+1)進行counts數轉換，接下來進行還原
a <- a^2-1
#a存儲的全是數值型向量，可以直接通過數學處理進行全表修改。
對表達矩陣進行篩選構建
實驗設計：列出TNBC的ID（前半部分），所有病人ID，有雙樣品類型的seq的病人ID，取TNBC的ID和雙樣品類型的seq病人ID交集作為目的病人ID。最後用 all_p %in% need_p 構建logical向量，被原始表達矩陣引用，構建符合目的要求的表達矩陣
##接下來是取樣本，需要取118個三陰性乳腺癌的樣本，並且是成對的樣本，即既有癌組織又有癌旁組織。
tnbc_p <- substring(tnbc_sample,1,12)
#取tnbc的每個元素的第一個到第12個字元
#tnbc是TNBC患者的ID前部分
all_p <- substring(colnames(a),1,12)
head(all_p)
#取表達矩陣的病理號前12位
table(all_p)
#這裡如果病人號前12位相同，說明是即既有癌組織又有癌旁組織。
table(all_p) == 2
paired_p <- names(table(all_p)[table(all_p) == 2])
#去除所有有兩個組織類型的病人ID前部分賦值給paired_id
need_p <- intersect(tnbc_p,paired_p)
#need_p是118個是成對的樣本。即既有癌組織又有癌旁組織
exprSet <- a[,all_p %in% need_p]
#構建表達矩陣：三陰性乳腺癌的樣本，並且是成對的樣本，即既有癌組織又有癌旁組織
#這裡通過 all_p %in% need_p 處理返回一個 logical向量，logical和第一個數據元素all_p一一對應，元素個數和順序和all_P相同。向量可以在在data.frame中引用，只會讀取 TRUE
對構建的表達矩陣進行數據篩選
實驗設計：利用apply函數統計條件符合情況，然後保存為logical向量，然後利用logical的被引用來篩選
table(apply(exprSet, 1, function(x){sum(x==0)<10}))
tmp <- apply(exprSet, 1, function(x){sum(x==0)<10})
#對exprSet進行按行執行函數function：
#如果表達矩陣是0則取1，在18個樣本中，如果這個基因的表達矩陣超過10個樣本都是空值，則捨棄這個基因
#tmp是一個logical向量
exprSet <- exprSet[tmp,]
save(exprSet,file = tnbc_paired_exprSet.Rdata)
#保存三陰性乳腺癌的雙樣本表達矩陣
癌症組織和正常組織的區分和標記
實驗設計：利用substr()函數截取患者ID的有效信息位置，用ifelse()函數進行判斷，並用factor()函數生成因子型變數。最後賦值得到患者的樣品類型
Tips：TCGA可以根據第14和第15位判斷是癌組織還是癌旁組織。01表示癌症組織，11表示正常組織
#TCGA可以根據第14和第15位判斷是癌組織還是癌旁組織。01表示癌症組織，11表示正常組織
colnames(exprSet)
#提取列名
substr(colnames(exprSet),14,15)
#提取列名字元串的第14和15位置字元
as.numeric(substr(colnames(exprSet),14,15))
#現在提取的字元變成數字 11和 11不等同
ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,tumor,normal)
#判斷第14和15號位的數值<10 即等於01時候，返回tumor 否則返回normal
factor(ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,tumor,normal))
#question 把它變成因子型，因子型是有順序的。這裡需要變換成因子型
group_list=factor(ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,tumor,normal))
#group_list記錄了表達矩陣中患者ID的對應的組織類型，和表達矩陣的順序一樣。
table(group_list)
表達矩陣的篩選（應該在上一步一起進行）
實驗設計：修正小於0的（因為變換過程中可能會產生-1，會影響實驗），檢查並去除na值。
exprSet <- ceiling(exprSet)
#?ceiling（）取整，不小於變數本身
#round()以四捨五入形式保存指定小數位數
#floor()不大於變數本身最大整數
#table(is.na(exprSet))
#pmax(),使得矩陣最小值是0？
save <- exprSet
exprSet[exprSet<0] <- 0
table(exprSet<0)
table(is.na(exprSet))
患者的癌症組織和正常組織的基因差異表達分析
library(DESeq2)
#構建一個病例號和腫瘤分類的對應關係
colData <- data.frame(row.names = colnames(exprSet),group_list= group_list)
#colData:是患者ID號和組織類型的對應關係
#構建DESeq()函數要求的表達式
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)
#countData = exprSet,指出DESeq的表達矩陣
#colData = colData,知名每個表達矩陣的分類，比如實驗組&對照組，正常組織&癌症組織
#question
#design= ~group_list，因子型，指出不同組的區別，是有順序的。
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~group_list)
dds <- DESeq(dds)
#進行差異基因分析
resultsNames(dds)
res <- results(dds, contrast=c("group_list","tumor","normal"))
#用group_list來做引導文件，用tumor來比較normal組織
resOrdered <- res[order(res$padj),]
#把res差異分析文件通過padj來排序
head(resOrdered)
resOrdered=as.data.frame(resOrdered)
#把resOrdered變成數據框，
write.table(resOrdered, file="single cell DGE order.xls", sep=" ",quote=F)
#寫出差異基因文件
基因的注釋
正在努力更新...
如果您感覺本文對您有些許幫助, 是筆者的莫大榮幸 .
推薦閱讀：
