TCGA_DESeq2分析

本篇是对“测序数据下载”、“DESeq2分析”的综合应用，因此会少一些说明：

“测序数据下载”详见 TCGA_改版后STAR-count处理方法_老实人谢耳朵的博客

“DESeq2分析”详见 TCGA_DESeq2分析_Gleason H vs L_老实人谢耳朵的博客

一、数据下载

library("TCGAbiolinks")
library("rjson")

1 下载tsv文件和json文件

tsv文件下载：

query <- GDCquery(project = "TCGA-PRAD", #项目名data.category = "Transcriptome Profiling",data.type = "Gene Expression Quantification",workflow.type = "STAR - Counts")
t_samplesDown <- getResults(query,cols=c("cases"))
#从sampleDown检索出TP和NT
t_dataSmTP <- TCGAquery_SampleTypes(barcode = t_samplesDown,typesample = "TP") #可选
t_dataSmNT =TCGAquery_SampleTypes(barcode = t_samplesDown,typesample = "NT") #可选
t_dataSmAll=c(t_dataSmTP,t_dataSmNT)
#筛选完成的query
queryDown <- GDCquery(project = "TCGA-PRAD",data.category = "Transcriptome Profiling",data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts", barcode = t_dataSmAll)
#下载GDCquery的结果
GDCdownload(queryDown,method = "api",directory = "GDCdata_star_count_TP&NT",files.per.chunk = 6) #网速慢可以设小一点

json文件下载：

详见 TCGA_改版后STAR-count处理方法_老实人谢耳朵的博客-CSDN博客

2 整合tsv文件、json文件到一个文件夹

3 提取count矩阵

#整理“all”中的文件
result <- fromJSON(file = "E:/R/PRAD Data Mining/PRAD_data_mining/TCGA/GDCdata_star_count_TP&NT/metadata.cart.2022-05-01.json")
metadata <- data.frame(t(sapply(result,function(x){id <-  x$associated_entities[[1]]$entity_submitter_idfile_name <- x$file_nameall <- cbind(id,file_name)
})))
rownames(metadata) <- metadata[,2]#获取example
example <- data.table::fread('E:/R/PRAD Data Mining/PRAD_data_mining/TCGA/GDCdata_star_count_TP&NT/all/005d2b9e-722c-40bd-aa5c-bd4e8842cb04.rna_seq.augmented_star_gene_counts.tsv',data.table = F)#读入一个tsv文件，查看需要的列数，“unstranded”#获取raw
t_dir <- 'E:/R/PRAD Data Mining/PRAD_data_mining/TCGA/GDCdata_star_count_TP&NT/all/'
t_samples=list.files(t_dir)
sampledir <- paste0(t_dir,t_samples) #各个文件路径raw <- do.call(cbind,lapply(sampledir, function(x){rt <- data.table::fread(x,data.table = F) #data.table::fread函数rownames(rt) <- rt[,1]rt <- rt[,4]###第4列为“unstranded”
}))#替换行名、列名
colnames(raw)=sapply(strsplit(sampledir,'/'),'[',8)###列名,'8'为文件名005d2b9e-722c-40bd-aa5c-bd4e8842cb04.rna_seq.augmented_star_gene_counts.tsv
rownames(raw) <- example$gene_id ##行名
raw_t <- t(raw)t_same <- intersect(row.names(metadata),row.names(raw_t)) dataPrep2 <- cbind(metadata[t_same,],raw_t[t_same,])
rownames(dataPrep2) <- dataPrep2[,1]
dataPrep2 <- t(dataPrep2)
dataPrep2 <-dataPrep2[-c(1:6),] #dataPrep2为未注释count矩阵

4 dataPrep2中数据类型为“character”，需要转为“numeric”

puried_data=apply(dataPrep2,2,as.numeric)

5 基因注释、预处理

#基因注释
rownames(puried_data)=example[5:nrow(example),'gene_name']#去除重复基因名
t_index=order(rowMeans(puried_data),decreasing = T)#计算所有行平均值，按降序排列
t_data_order=puried_data[t_index,]#调整表达谱的基因顺序
keep=!duplicated(rownames(t_data_order))#对于有重复的基因，保留第一次出现的那个，即行平均值大的那个
puried_data=t_data_order[keep,]#得到最后处理之后的表达谱矩阵#预处理
dataFilt =puried_data[rowMeans(puried_data)>10,] #剔除低表达基因
write.csv(dataFilt,file = "dataFilt__TP&NT.csv",quote = FALSE)

View(dataFilt)

二、DESeq2分析

rm(list=ls())library("DESeq2")

1 读入count矩阵——dataFilt

dataFilt=read.csv(file = "dataFilt_TP&NT.csv", header=T, row.names=1,check.names=FALSE)
dataFilt=as.matrix(dataFilt)

2 创建group分组矩阵——design

query <- GDCquery(project = "TCGA-PRAD", #项目名data.category = "Transcriptome Profiling",data.type = "Gene Expression Quantification",workflow.type = "STAR - Counts")
t_samplesDown <- getResults(query,cols=c("cases"))
t_dataSmTP <- TCGAquery_SampleTypes(barcode = t_samplesDown,                                   typesample = "TP") #可选
t_dataSmNT =TCGAquery_SampleTypes(barcode = t_samplesDown,typesample = "NT") #可选design=data.frame(c(t_dataSmTP,t_dataSmNT),c(rep('TP',length(t_dataSmTP)),rep('NT',length(t_dataSmNT))))
colnames(design)=c('Sample','Group')

View(design)

3 dds、dds1、res三步走

dds <- DESeqDataSetFromMatrix(countData = dataFilt, colData = design, design= ~Group) #countData需要一个count matrix，colData需要一个有分组信息的dataframe，design需要指定分组信息中的列
dds1 <- DESeq(dds,fitType = 'mean')
res <- results(dds1, contrast = c('Group', 'TP', 'NT')) #Group对应上面的“design= ~Group”，且TP在前，NT在后，表示TP比NT

4 获取差异比较结果

result <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)

View(result)

TCGA_DESeq2分析_TP vs NT相关推荐

TCGA_DESeq2分析_Gleason H vs L
一文掌握R包DESeq2的差异基因分析过程 - 简书 DESeq2包分析需两个矩阵: 1 count矩阵(必须为matrix) 行为Sample id 2 group矩阵(matrix或datafra ...
带有RESTEasy + JAXB + Jettison的JSON示例
RESTEasy使用Jettison JSON库在JSON之间来回映射JAXB注释对象. 在本教程中,我们向您展示如何将带有JAXB注释的对象转换为JSON格式并将其返回给客户端. 杰克逊(Jac ...
基因家族分析（3）进化树构建及美化
基因家族树构建最常用的方法是 NJ 法和 ML 方法,构建进化树之前,需要进行多序列比对. 欢迎关注Bioinfor 生信云微信公众号! 多序列比对多序列比对的作用是将核酸或者氨基酸序列对齐,使得相 ...
Wikidata知识图谱介绍与数据处理
1. Wikidata简介 Wikidata(维基数据)是一个自由开放的知识库,可以同时被人和机器阅读.编辑[1].根据官网介绍,Wikidata作为一种结构化数据的集中存储,为其他维基媒体(Wiki ...
什么是MBTI，16种人格类型详解
MBTI是Myers-Briggs Type Indicator的缩写,是一种经典的人格类型测评工具,通过对个体心理偏好的测试和分析. MBTI人格理论认为,人类的人格特征可以从四个角度进行分析:认知 ...
中小企业网管管理完全篇 [转]
现在,若有哪家企业还没有搭建起自己的网络系统,那么人们一定会说:"这家企业太落后了!"这种"落后"现象,虽然目前在我国的现代企业中不多见了,但是在占我国企业总数 ...
Oracle安装的一些问题收集[转]
在安装过程中出现的一些问题的解决办法.值得收藏与学习.比如在安装的时候如果有中文的路径则会出现类似这样的提示:加载数据库时出错:areasQueries Oracle的系统要求企业版:CPU 最低P ...
轮廓仪 wyko matlab,VEECO光学轮廓仪
VEECO光学轮廓仪 Wyko NT9100系列采用了Veeco第九代光学轮廓探测技术,在0.1nm到10mm的垂直扫描范围内提供了快速.高精度的三维表面形貌测量功能. NT9100光学轮廓仪是一款使 ...
c语言约瑟夫环分函数,c语言实现约瑟夫环问题
<c语言实现约瑟夫环问题>由会员分享,可在线阅读,更多相关<c语言实现约瑟夫环问题(16页珍藏版)>请在人人文库网上搜索. 1.一)基本问题1问题描述设有编号为1,2,小的n ...

TCGA_DESeq2分析_TP vs NT

一、数据下载

二、DESeq2分析

TCGA_DESeq2分析_TP vs NT相关推荐

最新文章

热门文章