實驗設計

實驗目的決定試驗方法和途徑拳亿。
試驗目的 :獲取三陰性乳腺癌的正常組織和癌癥組織的基因表達差異情況晴股，比較三陰性乳腺癌中的基因表達變化情況。
試驗設計 :通過TCGA獲取乳腺癌的RNA-seq表達數(shù)據(jù)肺魁，篩選出三陰性乳腺癌的樣本电湘，通過比較癌癥和正常組織的表達差異。

TCGA數(shù)據(jù)庫簡介

一句話介紹:TCGA數(shù)據(jù)庫是一個由國家癌癥研究所(National Cancer Institute)和美國人類基因組研究所(National Human Genome Research Institute)共同監(jiān)督的一個項目。使用對患者樣本的高通量基因組測序和分析技術來試圖提供括基因表達譜寂呛，拷貝數(shù)變異分析怎诫，SNP基因分型，全基因組DNA甲基化分析贷痪，微RNA分析等信息幻妓。收錄了33種癌癥基因組測序數(shù)據(jù)。TCGA數(shù)據(jù)處理和整理比Oncoman和GEO困難一些劫拢。但是針對腫瘤和癌癥所能提供的信息是很完善和可靠的肉津。
TCGA和GEO存在的區(qū)別是，GEO存在各種研究領域和研究方向的NGC數(shù)據(jù)和分析舱沧。TCGA是專門針對腫瘤和癌癥設立的妹沙。TCGA優(yōu)勢是豐富且規(guī)范的臨床數(shù)據(jù)，以及針對每種癌型的大樣本量熟吏。

TCGA數(shù)據(jù)的獲取

背景介紹：
三陰性乳腺癌是指癌組織免疫組織化學檢查結(jié)果為雌激素受體（ER）初烘、孕激素受體（PR）和原癌基因Her-2均為陰性的乳腺癌。這類乳腺癌占所有乳腺癌病理類型的10.0%～20.8%分俯，具有特殊的生物學行為和臨床病理特征肾筐，預后較其他類型差。--from:百度百科:三陰性乳腺癌

實驗設計：
三陰性乳腺癌的篩選標準是根據(jù)pheotype來確定的缸剪，表達量是通過RNAseq結(jié)果確定的吗铐，正常組織和癌癥組織是通過病例號確定的。

1. TCGA項目的數(shù)據(jù)可以通過Genomic Data Commons Data Portal獲取杏节，即通過GDC來訪問唬渗，訪問地址：https://portal.gdc.cancer.gov/。
2. TCGA數(shù)據(jù)庫公開免費奋渔，所以有許多針對TCGA數(shù)據(jù)進行整合的網(wǎng)站镊逝。還可以通過UCSC Xena進行下載：https://xena.ucsc.edu/public。試驗數(shù)據(jù)的選擇和目的息息相關嫉鲸，試驗設計：
   
   GDC
   
   點擊BRCA（乳腺癌）進入數(shù)據(jù)選擇界面撑蒜。
3. 選擇gene expression RNAseq>HTSeq - Counts。注意不要用RPKM等經(jīng)過了normlization的表達矩陣來分析玄渗。要使用Counts來進行差異分析座菠，因為在差異分析時候會自動進行標準化。如果數(shù)據(jù)經(jīng)過處理例如log2+1藤树，則可以下載后逆運算轉(zhuǎn)變回來浴滴。
4. 選擇phenotype>Phenotype這里面有病人的病例信息等，可以通過統(tǒng)計篩選出三陰性乳腺癌的患者岁钓。
   
   BRCA

5. 點擊連接進去會看到詳細的信息如下載地址升略、樣品數(shù)微王、數(shù)據(jù)處理方法等。

   
   
   RNAseq界面詳細信息

數(shù)據(jù)預處理

實驗設計：篩選出三陰性乳腺癌的患者ID品嚣，再篩選出同時有癌癥組織樣本和癌旁組織樣本骂远，計算初始Count值。

三陰性乳腺癌患者（TNBC）篩選

實驗設計：在R語言中處理數(shù)據(jù)腰根，選擇breast_carcinoma_estrogen_receptor_status（ER）激才、PR、HER2受體一欄全為隱形（Negative）的患者额嘿。
原始的phenotype文件如下圖瘸恼，信息量巨大

phenotype文件

讀取文件

phenotype_file <- read.table("TCGA-BRCA.GDC_phenotype.tsv",header = T, sep = '\t', quote = "")
phenotype_colnames <- as.data.frame(colnames(phenotype_file))
#讀取文件并去讀列名（penotype的類型）
table(phenotype_file$breast_carcinoma_estrogen_receptor_status)
#查看雌激素受體(ER)的表達狀態(tài)。
table(phenotype_file$breast_carcinoma_progesterone_receptor_status)
#查看孕激素受體（PR）狀態(tài)
table(phenotype_file$lab_proc_her2_neu_immunohistochemistry_receptor_status)
#查看原癌基因Her-2狀態(tài)

talbe

列出三項指標的列表册养，方便篩選

實驗設計：查找三項指標的列并列出

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號）

eph查看

函數(shù)學習：

?apply
# b為：
# first second
# one       1      4
# two       2      5
# three     3      6
# apply(b,1,sum)
# 上面的指令代表對矩陣b進行行計算靠闭，分別對每一行進行求和。函數(shù)涉及了三個參數(shù)：
# 第一個參數(shù)是指要參與計算的矩陣坎炼；
# 第二個參數(shù)是指按行計算還是按列計算愧膀，1——表示按行計算，2——按列計算谣光；
# 第三個參數(shù)是指具體的運算參數(shù)檩淋。
# 上述指令的返回結(jié)果為：
# one   two three 
# 5     7     9 

TNBC的篩選

實驗設計：查找指標的狀態(tài)并統(tǒng)計；列出三陰性并篩選出來這些數(shù)據(jù)

tnbc_rownum <- apply(eph,1,function(x)sum(x=="Negative"))
#eph：記錄三項指標類型的矩陣
#將eph按照列執(zhí)行sum(x=="Negative"):統(tǒng)計各行中萄金，eph列中有Negative的數(shù)量
#通過查看陰性的數(shù)量祝沸，函數(shù)返回結(jié)果是一個數(shù)字向量敛苇。內(nèi)容是0或1、2盖灸、3跛蛋。注意數(shù)字的排列順序是對應病人ID的行展懈。
tnbc_rownum
tnbc_sample <- phenotype_file[tnbc_rownum == 3, 1]
#通過讓 tnbc_rownum==3 判別橄浓，如果成立吠谢，會返回一個logical向量，logical向量是可以直接被矩陣引用的的圆，只會讀取TURE的行或著列
#統(tǒng)計每一行中數(shù)值是3的（即全為“Negative”）并取第一列即病人ID
#tnbc_sample：存儲著三陰性乳腺癌的
tnbc_sample
save(tnbc_sample,file = 'tnbc_sample.Rdata')
#保存Rdata

陰性指標統(tǒng)計

基因表達矩陣的構(gòu)建

基因表達矩陣的讀取和讀取后格式修改

實驗設計：轉(zhuǎn)換data.frame 鼓拧；行名修改半火；Count值還原

library(data.table)
#?data.table
#R中的data.table包提供了一個data.frame的高級版本越妈，讓你的程序做數(shù)據(jù)整型的運算速度大大的增加。
a <- fread("TCGA-BRCA.htseq_counts.tsv.gz",sep = '\t',header = T)
#讀取表達矩陣文件
a <- as.data.frame(a)
#轉(zhuǎn)換為數(shù)據(jù)框钮糖，讀取后梅掠，轉(zhuǎn)換前的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數(shù)阎抒，網(wǎng)站使用了log2(count+1)進行counts數(shù)轉(zhuǎn)換酪我，接下來進行還原
a <- a^2-1
#a存儲的全是數(shù)值型向量，可以直接通過數(shù)學處理進行全表修改且叁。

RNAseq結(jié)果的讀取與處理

對表達矩陣進行篩選構(gòu)建

實驗設計：列出TNBC的ID（前半部分）都哭，所有病人ID，有雙樣品類型的seq的病人ID逞带，取TNBC的ID和雙樣品類型的seq病人ID交集作為目的病人ID欺矫。最后用 all_p %in% need_p 構(gòu)建logical向量，被原始表達矩陣引用展氓，構(gòu)建符合目的要求的表達矩陣

##接下來是取樣本穆趴，需要取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]
#構(gòu)建表達矩陣：三陰性乳腺癌的樣本，并且是成對的樣本歪赢，即既有癌組織又有癌旁組織
#這里通過 all_p %in% need_p 處理返回一個 logical向量酪耕，logical和第一個數(shù)據(jù)元素all_p一一對應，元素個數(shù)和順序和all_P相同轨淌。向量可以在在data.frame中引用迂烁，只會讀取 TRUE

目的表達矩陣

對構(gòu)建的表達矩陣進行數(shù)據(jù)篩選

實驗設計：利用apply函數(shù)統(tǒng)計條件符合情況，然后保存為logical向量递鹉，然后利用logical的被引用來篩選

table(apply(exprSet, 1, function(x){sum(x==0)<10}))
tmp <- apply(exprSet, 1, function(x){sum(x==0)<10})
#對exprSet進行按行執(zhí)行函數(shù)function：
#如果表達矩陣是0則取1盟步，在18個樣本中，如果這個基因的表達矩陣超過10個樣本都是空值躏结，則舍棄這個基因
#tmp是一個logical向量
exprSet <- exprSet[tmp,]
save(exprSet,file = 'tnbc_paired_exprSet.Rdata')
#保存三陰性乳腺癌的雙樣本表達矩陣

癌癥組織和正常組織的區(qū)分和標記

實驗設計：利用substr()函數(shù)截取患者ID的有效信息位置却盘，用ifelse()函數(shù)進行判斷，并用factor()函數(shù)生成因子型變量媳拴。最后賦值得到患者的樣品類型

Tips：TCGA可以根據(jù)第14和第15位判斷是癌組織還是癌旁組織黄橘。01表示癌癥組織，11表示正常組織

#TCGA可以根據(jù)第14和第15位判斷是癌組織還是癌旁組織屈溉。01表示癌癥組織塞关，11表示正常組織
colnames(exprSet)
#提取列名
substr(colnames(exprSet),14,15)
#提取列名字符串的第14和15位置字符
as.numeric(substr(colnames(exprSet),14,15))
#現(xiàn)在提取的字符變成數(shù)字 '11'和 11不等同
ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,'tumor','normal')
#判斷第14和15號位的數(shù)值<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的（因為變換過程中可能會產(chǎn)生-1椰于，會影響實驗）怠益，檢查并去除na值。

exprSet <- ceiling(exprSet)
#?ceiling（）取整瘾婿，不小于變量本身
#round()以四舍五入形式保存指定小數(shù)位數(shù)
#floor()不大于變量本身最大整數(shù)
#table(is.na(exprSet))
#pmax(),使得矩陣最小值是0蜻牢？
save <- exprSet
exprSet[exprSet<0] <- 0
table(exprSet<0)
table(is.na(exprSet))

group_list和colData

患者的癌癥組織和正常組織的基因差異表達分析

library(DESeq2)
#構(gòu)建一個病例號和腫瘤分類的對應關系
colData <- data.frame(row.names = colnames(exprSet),group_list= group_list)
#colData:是患者ID號和組織類型的對應關系

#構(gòu)建DESeq()函數(shù)要求的表達式
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)
#countData = exprSet,指出DESeq的表達矩陣
#colData = colData,知名每個表達矩陣的分類，比如實驗組&對照組偏陪，正常組織&癌癥組織
#question
#design= ~group_list孩饼，因子型，指出不同組的區(qū)別竹挡，是有順序的镀娶。
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變成數(shù)據(jù)框揪罕，
write.table(resOrdered, file="single cell DGE order.xls", sep="\t",quote=F)
#寫出差異基因文件

DEG文件

基因的注釋

正在努力更新...