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、设置工作路径
setwd("C:\\Users\\86188\\Desktop\\CRC-TCGA\\UCSC xena\\1. UCSC Xena下载的数据")
3、导入表达谱数据
COADdata=read.table("TCGA-COAD.htseq_counts.tsv",header=T,sep='\t')
COADdata[1:4,1:4]
4、去除Ensembl_ID小数点后面的数字
COADdata1<-separate(COADdata,Ensembl_ID,into= c("Ensembl_ID"),sep="\\.")
COADdata1[1:4,1:4]
5、下载Homo_sapiens.GRCh38.99.chr.gtf.gz文件，然后放到我们R语言的工作目录内，打开文件（可以在网上搜到下载方法，此不在详述），进入相应的工作路径
setwd("C:\\Users\\86188\\Desktop\\CRC-TCGA\\UCSC xena\\2. Ensemble_ID转换为Gene_Symbol")
gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz')
class(gtf)
6、转换为数据框
gtf <-as.data.frame(gtf)
7、查看文件，保存文件为Rdata，将来方便我们直接打开
dim(gtf)
save(gtf,file ="Homo_sapiens.GRCh38.99基因组注释文件.Rda")
8、#由于GTF文件中，基因ID的列名是gene_id,所以我们把COADdata1中的基因列名改成一样的，方便后续合并
colnames(COADdata1)[1]<-'gene_id'
9、通过浏览文件看到我们需要的主要信息在两列：1. type，这一列我们需要选择gene，2. gene_biotype，这一列我们需要选择protein_coding，当然你也可以选择其他的种类，比如miRNA，长链非编码等。我们首先把蛋白编码的基因的行都筛选出来
a=dplyr::filter(gtf,type=="gene",gene_biotype=="protein_coding")
dim(a)
10、这个时候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函数）
d=select(c,-gene_id,-gene_biotype)，mRNAdata=distinct(d,gene_name,.keep_all = T)
dim(mRNAdata)
13、把行名由数字换成基因
rownames(mRNAdata)<-mRNAdata[,1]
mRNAdata<-mRNAdata[,-1]
14、我们下载的数据取过了log2（count+1），这里我们再返回count
mRNAdata<-2^mRNAdata -1
View(mRNAdata)
class(mRNAdata)#此时是数据框“data.frame”，与“matrix”区别就是前者行名不允许有重复的
15、#如果某基因在样本间的整体表达水平很低，则小幅度的read counts扰动即可造成较大的fold change差异，但这种差异属于noise，并不具有生物学意义。所以需要设置一定的threshold将这些基因去除，减少假阳性。过滤的原则：原理类似，具体操作各自不同。本例使用相对更苛刻的原则：若某基因在大于80%的样本中read counts数<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”，即空格，制表符，换行符 。到此为止，差异分析前的基因表达数据就整理好了
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个列名为空
save(mRNAdata,file ="mRNAdata.Rda")
