R语言 2022 TCGA数据库转录组提取 新版TCGA 表格提取 一键精灵

一键生成版本 （找工作去了，不改了）

# 1. 下载数据与json文件

#  2. 不同文件夹文件提取
# 将次级目录的文件夹里面的文件提取到同一个文件夹下

# 一些基础操作
list.files(pattern = "\\.tsv") #当前目录显示以tsv结尾的文件
dir()
dir(all.files=TRUE)  #显示隐藏的文件
getwd()

# 验证二级文件夹里面是否有文件，无返回NA值，为TRUE，
filename=1
for (i in 1:200) { 
  a=as.character(list.files(list.files()[i])[1])
  ifelse( a%in% NA==TRUE, NA,'b')
  b=paste(getwd(),"/",list.files()[i],"/",list.files(list.files()[i])[1],sep = "")
  filename[[i]]=b
}
filename
filename <- as.data.frame(filename)
# 对200行数据进行过滤
filename2=apply(filename, 2, 
                function(x){gsub(pattern = ".*(NA).*", 
                                 replacement = "\\21",x) })
filename2 <- as.data.frame(filename2)
filename2 <- filename2[filename2$filename!=1,]
filename2 #这时候提取到试运行的10个文件的绝对路径

#复制文件到同一文件夹
#获得路径，在后面加入/data  
# dir.create('C:/Users/shao/Desktop/T2/data3') #创建一个目录
eval(parse(text = paste0("dir.create('",getwd(),"/data')")))
# 复制文件到data这个文件夹
file.copy(filename2,"data") 


# 先读取jason文本信息,此时工作路径还没转向/data
library(rjson)
result <- fromJSON(file="metadata.cart.2022-05-04.json") #先读取
# 转向data目录
eval(parse(text = paste0("setwd('",getwd(),"/data')")))
getwd()


# 3. 提取json的信息
# 提取json的信息，匹配文件名与对应的TCGA- - - - 情况
Metadata=data.frame()
for (i in 1:200) {
  a <- result[[i]][["file_name"]]
  b <- result[[i]][["associated_entities"]][[1]][["case_id"]]
  c <- result[[i]][["associated_entities"]][[1]][["entity_submitter_id"]]
  Metadata[i,1] =a
  Metadata[i,2] =b
  Metadata[i,3] =c
}
names(Metadata)
names(Metadata)[1] <- "file_name"
names(Metadata)[2] <- "case_id"
names(Metadata)[3] <- "entity_submitter_id"
names(Metadata)
# 查看重复情况 可以预推断重复情况
# PS:# 有的病例没有对照/只有对照情况
table(duplicated(Metadata$case_id))
write.csv(Metadata,"metadataID1.csv")


# 4. 批量改文件名

# 都是在EXCEL生成相应的代码
# 4.1 *linux改* Git bash here 修改路径的文件名
# mv file_name.tsv TCGA- - - .tsv
# 4.2 或者根据生成的metadataID1.csv，修改R里面的data
# TCGA_CN_5360_01A_01R_1436_07.tsv <- 6303629f-y.tsv


# 4.3 读取的时候就修改

# 按新名字读取数据

# 提取新文件名
# 读取数据
lf <-list.files(pattern = ".tsv$") #以report.tsv 结尾的
files <- gsub("\\.tsv", "", lf)   #切掉后缀.tsv，获得这些名称，为循环准备
files 
# 一定要一一对应，或可以利用metadataID1
metadataID1 <- read.csv("metadataID1.csv")
# 【先将 -改成_】
metadataID1$entity_submitter_id <- gsub(pattern = "\\-",
                   replacement = "\\_",metadataID1$entity_submitter_id ) 
metadataID1$entity_submitter_id[1:5]

lf <- as.data.frame(lf)
lf
metadataID2 <- dplyr::full_join(lf,metadataID1,by=c("lf"="file_name"))
nrow(metadataID2)
for (i in 1:123){
  assign(metadataID2$entity_submitter_id[i], 
         read.table(metadataID2$lf[i], sep = '\t', header = TRUE))
}

# 5. 循环读取.tsv [已做4.3跳过]
# lf <-list.files(pattern = ".tsv$") #以report.tsv 结尾的
# files <- gsub("\\.tsv", "", lf)   #切掉后缀.tsv，获得这些名称，为循环准备
# files
# for (i in seq_along(files))
#   assign(files[i], read.table(lf[i], sep = '\t', header = TRUE))


# 6. 提取要的部分
# TCGA_CN_5360_01A_01R_1436_07 = TCGA_CN_5360_01A_01R_1436_07[-c(1:4), c("gene_id","fpkm_unstranded")]
# names(TCGA_CN_5360_01A_01R_1436_07)[2] <- "TCGA_CN_5360_01A_01R_1436_07"


# 借助ECXEL提取
# 新生成代码,导出复制到新Untitled执行

# 提取要的数据
files2=ls(pattern = "TCGA")
# 删掉空白的1：4 提取要的两列
for (i in 1:10) {
  eval(parse(text = paste0(files2[i], "<-", 
            files2[i],"[-c(1:4),c('gene_id','fpkm_unstranded')]")))
}

# 改名字
files2=ls(pattern = "TCGA");files2

for (i in 1:10) {
  eval(parse(text = paste0("names(",files2[i], ")[2]", "=",
                          "'",files2[i],"'")))
}

# 7. 多个数据合并
multimerge<-function(dat=list(),...){
  if(length(dat)<2)return(as.data.frame(dat))
  mergedat<-dat[[1]]
  dat[[1]]<-NULL
  for(i in dat){
    mergedat<-merge(all=TRUE,mergedat,i,...)
  }
  return(mergedat)
}

files2=ls(pattern = "TCGA");files2

listALL=list()

# 数据框合成list
for (i in 1:10) {
  eval(parse(text = paste0("listALL","[[",i,"]]", " <- ",files2[i])))
}

dataALL <- multimerge(listALL)
# 导出基因矩阵
write.csv(dataALL,"dataALL.csv",row.names = F)



# 8. 基因注解过程

# [下载基因信息文件](http://asia.ensembl.org/index.html)
# [具体下载位置](http://asia.ensembl.org/Homo_sapiens/Info/Index)
#下载 Homo_sapiens.GRCh38.106.gtf.gz 

# 8.1 读取文件GTF文件
# BiocManager::install("rtracklayer")
#   library("rtracklayer")
#   gtf_data = import('Homo_sapiens.GRCh38.106.gtf') 
#   gtf_data = as.data.frame(gtf_data)
#   names(gtf_data)
#   a=head(gtf_data,100)
#   write.csv(a,"GTF预览.csv")

# 8.2 处理未注解的基因矩阵文件

# 保留.前面的数字 
# ENSG基因后面有带 点. 为版本号
suppressMessages(library(tidyverse))
table(duplicated(dataALL$gene_id))  # 唯一ID
dataqie=dataALL
dataqie <- dataqie   %>% 
  separate(gene_id,  sep = "\\.",    #切割点要加\\
           into = c("gene_id2", "deleteVar")) %>%  
  select(-deleteVar)
dataqie[1:5,]
table(duplicated(dataqie$gene_id2))


# 查看重复的情况
a=dataqie[duplicated(dataqie[ ,'gene_id2']), ]
a
name <- t(a[,'gene_id2'])
dd <- as.vector(name)
dd
b=dataqie[which(dataqie$gene_id2%in% dd),]
b[, c(1:7)]

# 可以计算下重复的有没有表达
# 没有直接去掉
apply(b[,2:11],1,sum) #不处理也可以,后续的基因过滤也会去掉



# 8.3 处理与提取注解文件
#    table(duplicated(gtf_data$gene_id)) #是重复的
#    table(gtf_data$gene_biotype)
#    #查看重复是不是有唯一的gene biotype
#    b=table(gtf_data$gene_id,gtf_data$gene_biotype)
#    b[b>0] <- 1
#    b2 <- as.data.frame(b)
#    
#    b2[b2$Freq>1,] #检测有无大于1的，相当于有重复
#    #<0 行> (或0-长度的row.names) 代表gene_id与gene_biotype完全一一对应
#    table(b2$Freq)
#    #保留Freq=1的
#    table(duplicated(b2$Var1))
#    gtf=b2[b2$Freq==1,]
#    dim(gtf)
#    names(gtf)[1] <- "gene_id"
#    names(gtf)[2] <- "gene_biotype"
#    gtf <- gtf[,-3]
#    str(gtf)
#    gtf$gene_id <- as.character(gtf$gene_id)
#    str(dataqie$gene_id2)
#    # 提取注解文件
#    # 先查看gene_id和gene_name是不是一一对应
#    gtf_data$idandName <- paste(gtf_data$gene_id,"_",gtf_data$gene_name)
#    table(duplicated(gtf_data$idandName))
#    nrow(gtf) #与上面唯一值对应，证明可以直接合并
#    
#    #
#    #检查ID重复情况
#    table(duplicated(gtf_data$gene_name))
#    #删除重复的行
#    gtf_data2 <- gtf_data[!duplicated(gtf_data$gene_name), ]
#    #在检查下
#    table(duplicated(gtf_data2$gene_name))
#    
#    
#    # 合并
#    gtf3 <- dplyr::left_join(gtf,gtf_data2,by=c('gene_id'='gene_id'))
#    is.na(gtf3$gene_name)
#    names(gtf3)
#    # 挑选要的变量
#    
#    gtffinal <- gtf3[ ,c('gene_id','gene_biotype.x',"gene_name")]
#    write.csv(gtffinal,"gtf最终注解文件")

#write.table(lncRNA.gtf,file = "lncRNA.gtf",sep = "\t",col.names = T)

# 这时候基因表达矩阵文件为 dataqie - dataqie$gene_id2
# 注解文件为gtffinal -

gtffinal <- read.csv("gtf最终注解文件.csv")
# 基因注解，也就是合并文件
GeneMatrix <- dplyr::left_join(dataqie,gtffinal,by=c('gene_id2'='gene_id'))