1、安裝并加載部分R包
#if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
#if(!require("rtracklayer")) BiocManager::install("rtracklayer")
#if(!require("tidyr")) BiocManager::install("tidyr")
#if(!require("dplyr")) BiocManager::install("dplyr")
#if(!require("DESeq2")) BiocManager::install("DESeq2")
#if(!require("limma")) BiocManager::install("limma")
#if(!require("edgeR")) BiocManager::install("edgeR")
library("rtracklayer")
library("tidyr")
library("dplyr")
2、設(shè)置工作路徑
setwd("C:\\Users\\86188\\Desktop\\CRC-TCGA\\UCSC xena\\1. UCSC Xena下載的數(shù)據(jù)")
3、導(dǎo)入表達(dá)譜數(shù)據(jù)
COADdata=read.table("TCGA-COAD.htseq_counts.tsv",header=T,sep='\t')
COADdata[1:4,1:4]
4、去除Ensembl_ID小數(shù)點(diǎn)后面的數(shù)字
COADdata1<-separate(COADdata,Ensembl_ID,into= c("Ensembl_ID"),sep="\\.")
COADdata1[1:4,1:4]
5、下載Homo_sapiens.GRCh38.99.chr.gtf.gz文件，然后放到我們R語言的工作目錄內(nèi)，打開文件（可以在網(wǎng)上搜到下載方法，此不在詳述），進(jìn)入相應(yīng)的工作路徑
setwd("C:\\Users\\86188\\Desktop\\CRC-TCGA\\UCSC xena\\2. Ensemble_ID轉(zhuǎn)換為Gene_Symbol")
gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz')
class(gtf)
6、轉(zhuǎn)換為數(shù)據(jù)框
gtf <-as.data.frame(gtf)
7、查看文件，保存文件為Rdata，將來方便我們直接打開
dim(gtf)
save(gtf,file ="Homo_sapiens.GRCh38.99基因組注釋文件.Rda")
8、#由于GTF文件中，基因ID的列名是gene_id,所以我們把COADdata1中的基因列名改成一樣的，方便后續(xù)合并
colnames(COADdata1)[1]<-'gene_id'
9、通過瀏覽文件看到我們需要的主要信息在兩列：1. type，這一列我們需要選擇gene，2. gene_biotype，這一列我們需要選擇protein_coding，當(dāng)然你也可以選擇其他的種類，比如miRNA，長鏈非編碼等。我們首先把蛋白編碼的基因的行都篩選出來
a=dplyr::filter(gtf,type=="gene",gene_biotype=="protein_coding")
dim(a)
10、這個(gè)時(shí)候a文件只有19939行了，下來我們只選擇gene_name，gene_id和gene_biotype這三列，其他都不要了
b=dplyr::select(a,c(gene_name,gene_id,gene_biotype))
b[1:4,]
11、接下來用COADdata1和b文件中共有的gene_id列合并文件
c=dplyr::inner_join(b,COADdata1,by ="gene_id")
c[1:5,1:5]
12、再去掉第2列(gene_id)，3列(gene_biotype)，基因名再去重（distinct函數(shù)）
d=select(c,-gene_id,-gene_biotype)，mRNAdata=distinct(d,gene_name,.keep_all = T)
dim(mRNAdata)
13、把行名由數(shù)字換成基因
rownames(mRNAdata)<-mRNAdata[,1]
mRNAdata<-mRNAdata[,-1]
14、我們下載的數(shù)據(jù)取過了log2（count+1），這里我們?cè)俜祷豤ount
mRNAdata<-2^mRNAdata -1
View(mRNAdata)
class(mRNAdata)#此時(shí)是數(shù)據(jù)框“data.frame”，與“matrix”區(qū)別就是前者行名不允許有重復(fù)的
15、#如果某基因在樣本間的整體表達(dá)水平很低，則小幅度的read counts擾動(dòng)即可造成較大的fold change差異，但這種差異屬于noise，并不具有生物學(xué)意義。所以需要設(shè)置一定的threshold將這些基因去除，減少假陽性。過濾的原則：原理類似，具體操作各自不同。本例使用相對(duì)更苛刻的原則：若某基因在大于80%的樣本中read counts數(shù)<10，則將其濾除。
qualified_genes <- c()
for (genes_in_sheet in rownames(mRNAdata)) {
qualification <- mRNAdata[genes_in_sheet,] <= 10
if (sum(qualification) < 0.8*length(mRNAdata)) {
qualified_genes <- append(qualified_genes,genes_in_sheet)
}
}
mRNAdata <- mRNAdata[qualified_genes,]
dim(mRNAdata)
16、保存文件，常見空白分隔符有：sep=” ”；sep = “\t”；sep = “\n”，即空格，制表符，換行符 。到此為止，差異分析前的基因表達(dá)數(shù)據(jù)就整理好了
write.csv(mRNAdata,"mRNAdata.csv",quote = F,row.names = T)
write.table(mRNAdata, file = "mRNAdata.txt", sep = "\t",row.names = T,col.names = NA,quote = F)#col.names = NA表示第1個(gè)列名為空
save(mRNAdata,file ="mRNAdata.Rda")
