m6A-lncRNA week01

title: "m6A-LGG-week1"
author: "Hannisal"
date: "2020/10/3"

knitr::opts_chunk$set(echo = TRUE)
options(Bioc_mirror="http://mirrors.ustc.edu.cn/bioc/")
options("repos"=c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if(!require(stringr))install.packages("stringr")
if(!require(rtracklayer))BiocManager::install("rtracklayer")
if(!require(stringr))install.packages("rlang")
if(!require(dplyr))install.packages("dplyr")
if(!require(tidyr))install.packages("tidyr")
library("tidyr") 
library("dplyr")
library("stringr")
library("rtracklayer")

1.TCGA的数据下载

用TCGA官网的GDC工具下载，需要准备gdc-client.exe和gdc_manifest文件

https://portal.gdc.cancer.gov/linked phrase

knitr::opts_chunk$set(eval = F)
options(stringsAsFactors = F)

cancer_typer="TCGA-LGG"
# 创建一个文件夹，用来储存下载的表达数据
if(!dir.exists("expdata"))dir.create("expdata")

# running in Terminal for windows
# ./gdc-client.exe download -m gdc_manifest_20200917_144518.text -d expdata

2.读取并整合数据

因为TCGA下载下来的表达数据，每个样本是一个压缩包一个txt文档，所以需要循环读取

knitr::opts_chunk$set(eval = F)
# 先任选两个counts文件读取，并观察geneid的顺序是否一致。
options(stringsAsFactors = F)
x = read.table("expdata/0aaaa7cc-b59d-4e08-ba7f-14e2d2ae4860/63a70793-af6d-4248-9085-5d69276a3fff.FPKM.txt.gz")
x2 = read.table("expdata/0b4334ce-bcab-4aed-a07b-0ccf6766d236/430accb2-dcaa-4d25-9405-e8ef97126eb4.FPKM.txt.gz")
identical(x$V1,x2$V1)
#顺序是一致的，可以直接cbind，不会导致顺序错乱。
#批量读取所有的counts.gz文件。
count_files = dir("expdata/",pattern = "*.FPKM.txt.gz$",recursive = T)

ex = function(x){
  result <- read.table(file.path("expdata/",x),row.names = 1,sep = "\t") 
  return(result)
  }
exp = lapply(count_files,ex)
exp <- do.call(cbind,exp)
dim(exp)
exp[1:4,1:4]
#发现问题：这样产生出来的表达矩阵没有列名。
#解决办法：找到一个文件名与样本ID一一对应的文件。cart-json文件。
meta <- jsonlite::fromJSON("metadata.cart.2020-09-18.json")
colnames(meta)
meta$associated_entities[1]

ids <- meta$associated_entities;class(ids)
ids[[1]]

class(ids[[1]][,1])

#可以看到，meta$associated_entities是个列表，这个列表里包含数据框，数据框的第一列内容就是tcga样本id。
#注意，换了数据需要自己探索存放在哪一列。不一定是完全一样的，需要确认清楚。
colnames(ids)
ID = sapply(ids,function(x){x[,4]})
file2id = data.frame(file_name = meta$file_name, ID = ID)

# 文件名与TCGA样本ID的对应关系已经得到，接下来是将其添加到表达矩阵中，成为行名。需要找到读取文件的顺序，一一对应修改。
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
count_files2[1] %in% file2id$file_name
file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
exp[1:4,1:4]

# save(exp,file = "exp.Rdata")

3.整理临床信息

# 1.根据样本ID整理临床信息
## 读入CiBioportal网站下载的TCGA的LGG临床数据
cl<- read.table("lgg_tcga_clinical_data.tsv",sep="\t",header = T,check.names = F)

ovtime<-data.frame(cl$`Patient ID`,cl$`Sample ID`,cl$`Overall Survival (Months)`,cl$`Overall Survival Status`)
colnames(ovtime)<-c("Patient_ID","Sample_ID","OStime","OSstatus")

knitr::kable(head(ovtime),caption = "ovtime")

## 查看样本重复情况
sum(duplicated(ovtime$Patient_ID))
sum(duplicated(ovtime$Sample_ID))

## 查看TCGA_ID的14，15位，去除转移样本，只保留原发灶,因为都来源于同一位患者
ovtime$Sample_ID<-substr(ovtime$Sample_ID,14,15)
table(ovtime$Sample_ID)
ovtime01<-ovtime[ovtime$Sample_ID=="01",]
## 按照文中的筛选条件：missing os values and os < 30days excuded
ovtime01<-ovtime01[!is.na(ovtime01$OStime),]
ovtime01<-ovtime01[!is.na(ovtime01$OSstatus),]
ovtime01<-ovtime01[!ovtime01$OStime<1,]

ovtime01$OSstatus[ovtime01$OSstatus=="0:LIVING"]=0
ovtime01$OSstatus[ovtime01$OSstatus=="1:DECEASED"]=1

knitr::kable(head(ovtime01),caption = "ovtime")

# save(ovtime01,file = "clinic.Rdata")



# 2.根据样本ID整理表达矩阵以及匹配临床信息和表达矩阵的样本ID
## 查看exp里的样本重复情况,因为有些病人可能收集了多个部位的，或者有转移样本
sum(duplicated(substr(colnames(exp),1,12)))
table(substr(colnames(exp),14,15))
# 可以看出这18个都是转移灶"02"
# 去除02的，再看看
exp<-exp[,substr(colnames(exp),14,15)=="01"]

## 把ENTRIZ ID去掉小数点以及TCGA_ID的后几位
colnames(exp)<-substr(colnames(exp),1,12)
rownames(exp)<-substr(rownames(exp),1,15)

## 匹配有临床信息的样本
sum(colnames(exp)%in%ovtime01$Patient_ID)
exp_match<-exp[,colnames(exp)%in%ovtime01$Patient_ID]



lnc_gtf<-rtracklayer::import("gencode.v35.long_noncoding_RNAs.gtf")
lnc_gtf<-as.data.frame(lnc_gtf)

head(lnc_gtf)

# save(exp_match,file = "exp_match.Rdata")
# save(ovtime01,file = "ovtime01.Rdata")


## 从lnc_gft提取需要的信息
lnc_id<- lnc_gtf[,10:12]
table(lnc_id$gene_type)
## TEC表示有待实验验证的片段

## 提取lncRNA对应的基因名
lncid<-lnc_id[lnc_id$gene_type=="lncRNA",] # 只保留lncRNA
lncid<-lnc_id # lncRNA和TEC
lncid$gene_id<-substr(lncid$gene_id,1,15)
lncid<-lncid[!duplicated(lncid$gene_id),]
sum(duplicated(lncid$gene_name)) # gene_name还是有9个重复
# 查看那9个鬼东西
gdx<-lncid[duplicated(lncid$gene_name),]
gdx1<-lncid[lncid$gene_name%in%gdx$gene_name,]

head(gdx1) ## 可以看到，有9个LNC是多个ensembl号对应一个基因名

lncid<-lncid[!duplicated(lncid$gene_name),]
## 按照lnc注释提取表达矩阵
lncid<-lncid[lncid$gene_id%in%rownames(exp_match),]
exp_match<-exp_match[rownames(exp_match)%in%lncid$gene_id,]


exp_match$gene_id<-rownames(exp_match)
exp_match$gene_name<-lncid$gene_name[match(rownames(exp_match),lncid$gene_id)]
## 先注释，再去掉重复的基因名
exp_match<-exp_match[!duplicated(exp_match$gene_name),]
rownames(exp_match)<-exp_match$gene_name
exp_match<-exp_match[,-((ncol(exp_match)-1):ncol(exp_match))]
## 删除都为0的基因
sum(apply(exp_match, 1, mean)==0)
exp_match<-exp_match[apply(exp_match, 1, mean)!=0,]

save(exp_match,file="exp_genename.Rdata")

问题：多个ensembl ID对应一个基因名，是取均值或者最大值？我们是随便取的其中一个。
结果：因为从官网下载的lncRNA注释文件里边有TEC和LNC的，TEC是指需要实验验证的片段，只用LNC提取了13604个，两者都用有14592

初步得到：

477个有临床信息和表达矩阵的LGG样本（表达矩阵保存为exp_match.Rdata），exp_genename.Rdata是提取的lncRNA表达矩阵，ovtime01是477个样本的生存信息。

4.提取21个m6a相关基因的表达矩阵

## 提取m6A相关基因的表达矩阵,需要重新导入exp_match，因为之前提取lncRNA改变了exp_match
m6A_id<-read.csv("m6a.csv",header = T,sep = ",",check.names = F)
exp_m6a<-exp_match[rownames(exp_match)%in%m6A_id$`Gene stable ID`,]
rownames(exp_m6a)<-m6A_id$`Gene name`[match(rownames(exp_m6a),m6A_id$`Gene stable ID`)]

save(exp_m6a,file = "exp_m6a.Rdata")

5.批量作相关性分析

准备两个文件，exp_m6a,exp_match(lncRNA)，一个是m6a的表达矩阵，一个是lncRNA的表达矩阵，用for循环对每个m6a基因和lncRNA进行相关性分析


## 准备两个文件，exp_m6a,exp_match(lncRNA)
exp_match<-t(exp_match)
exp_m6a<-t(exp_m6a)


## for循环
cor_result<-list()
for (i in 1:ncol(exp_m6a)){
  exp_cor<-cbind(exp_m6a[,i],exp_match)
  colnames(exp_cor)[1]<-colnames(exp_m6a)[i]
  
  ## exp_cor已就绪
  y <- as.numeric(exp_cor[,1])
  colnames <- colnames(exp_cor) 
  cor_data_df <- data.frame(colnames)
  for (m in 1:length(colnames)){
    test <- cor.test(as.numeric(exp_cor[,m]),y,type="pearson") 
    cor_data_df[m,2] <- test$estimate
    cor_data_df[m,3] <- test$p.value
  }
  names(cor_data_df) <- c("symbol","correlation","pvalue")
  ## pvalue<0.001&abs(correlation)>0.5
  cor_data_sig <- cor_data_df %>%
    filter(pvalue < 0.001) %>%
    filter(abs(correlation)>0.5)
  cor_result[[i]]<-cor_data_sig
  print(i)
}

save(cor_result,file = "cor_result.Rdata")
cor_result[[1]]

## 将相关性分析后得到的基因名提取出来
cor_gene<-matrix(data = NA,ncol = 1)
for (i in 1:21) {
  mt<-cor_result[[i]]$symbol[-1]
  mt<-as.matrix(mt[!mt%in%cor_gene])
  cor_gene<-rbind(cor_gene,mt)
  print(i)
}
cor_gene<-cor_gene[-1,]
length(unique(cor_gene))
## 716个
cor_gene_num<-unique(cor_gene)
## 保存txt文件
write.table(cor_gene_num,file = "cor_gene_num.txt")

结果

保存到一个list里面（cor_result.Rdata)，包含每个m6a基因对应的相关性分析结果（系数和P值），所有lncRNA取并集得到cor_gene_num.txt

6.单因素cox分析

477个TCGA的LGG样本，716个与m6a相关的lncRNA

## 数据准备：716个lncRNA的表达矩阵和之前匹配好的临床信息
lnc<-read.table("cor_gene_num.txt",header = T,sep = " ")

exp_lnc<-t(exp_match)
exp_lnc<-exp_lnc[,colnames(exp_lnc)%in%lnc$x]
## 提取477个有表达矩阵的样本的临床信息
ovtime01<-ovtime01[match(rownames(exp_lnc),ovtime01$Patient_ID),]
## 为单因素cox分析做准备
unicox<-cbind(ovtime01,exp_lnc)

## 创建一个空白的数据框用于存储分析结果
univar_out <- data.frame(matrix(NA,ncol(unicox)-4,6))
univar_out[,1]<-colnames(unicox)[-(1:4)]
colnames(univar_out)<-c("Genesymbol","Coeffcient","HR","lower .95","upper .95","P-value")

## cox_data即为输入文件,关于为什么不用unicox，因为matrix把数据类型给弄乱了而且相关性分析需要数据服从方差齐正太分布特点
cox_data<-unicox[,c(5:ncol(unicox))]
class(cox_data)
cox_data<-as.matrix(cox_data)
cox_data<-apply(cox_data, c(1,2),as.numeric)

# cox_data<-log2(cox_data + 0.01)# 左偏转正太分布，好像不用转换
cox_data<-cbind(unicox[,1:4],cox_data)

cox_data$OSstatus<-as.numeric(cox_data$OSstatus)

library(survival)
for (i in 1:(ncol(cox_data)-4)){
  cox = coxph(Surv(OStime,OSstatus) ~ cox_data[,i+4],data = cox_data)
  cox_summ = summary(cox)
  univar_out[i,2]=cox_summ$coefficients[,1]
  univar_out[i,3]=cox_summ$coefficients[,2]
  univar_out[i,4]=cox_summ$conf.int[,3]
  univar_out[i,5]=cox_summ$conf.int[,4]
  univar_out[i,6]=cox_summ$coefficients[,5]
}
save(univar_out,file = "univar_out.RData")
univar_out_0.05 = univar_out[univar_out[,6]<0.05,]
write.table(univar_out_0.05,file = "univar_out_0.05.txt",sep = "\t",quote = F)

单因素cox得到了313个lncRNA，保存在Univar_out_0.05.txt

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,271评论 5赞 476
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,275评论 2赞 380
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,151评论 0赞 336
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,550评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,553评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,559评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,924评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,580评论 0赞 257
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,826评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,578评论 2赞 320
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,661评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,363评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,940评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,926评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,156评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,872评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,391评论 2赞 342