二、UCSC Xena数据库下载表达数据Ensemble_ID转换为Gene并整理数据

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")
最后编辑于
?著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,172评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,346评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 159,788评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,299评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,409评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,467评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,476评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,262评论 0 269
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,699评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,994评论 2 328
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,167评论 1 343
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,827评论 4 337
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,499评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,149评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,387评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,028评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,055评论 2 352

推荐阅读更多精彩内容