一键生成版本 （找工作去了，不改了）
# 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:
# [具体下载位置](http:
#下载 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'))
扫一扫
290
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
