专栏首页生信小驿站单基因生信分析流程(6)单基因相似性分析

单基因生信分析流程(6)单基因相似性分析

第一步,下载COAD数据
##########################################################################################
## step1 load package and change Working Directory
###########################################################################################

library(TCGAbiolinks)
library(dplyr)
library(tidyr)
library(tibble)
library(edgeR)
library(limma)
rm(list=ls())
setwd('D:\\SCIwork\\F20ELFN1\\COAD')


##########################################################################################
## step2 download the expresssion data of lncRNA and mRNA
###########################################################################################


query <- GDCquery(project = "TCGA-COAD", 
                  data.category = "Transcriptome Profiling", 
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts")



GDCdownload(query, method = "api", files.per.chunk = 50)

library(SummarizedExperiment)

expdat <- GDCprepare(query = query,   save = TRUE, save.filename = "exp.rda")


count_matrix = as.data.frame(assay(expdat))
第二步,注释表达量
##########################################################################################

###########################################################################################



rm(list=ls())

setwd('D:\\SCIwork\\F20ELFN1\\COAD')

load('exp.rda')

count_matrix = as.data.frame(assay(data))

count_matrix[1:4,1:4]



fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}


expr <- as.data.frame (apply(count_matrix  , 2, fpkmToTpm))


expr <-  expr %>% rownames_to_column("gene_id")




##########################################################################################

###########################################################################################





setwd("D:\\Originaldata\\GRCH\\Homo_sapiens.GRCh38.90")

load("gtf_df.Rda")

test <- gtf_df[1:5,]

View(test)




mRNA_exprSet <- gtf_df %>%
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>%
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%
  dplyr::inner_join(expr,by ="gene_id") %>%
  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")


save(mRNA_exprSet,file = "mRNA_exprSet.Rda")


mRNA_exprSet <- mRNA_exprSet %>%
  tidyr::separate(gene_id, c("gene_name","gene_id","gene_biotype"),
                  sep = " \\| ")


mRNA_exprSet <- mRNA_exprSet[,-(2:3)]


index <- duplicated(mRNA_exprSet$gene_name)
mRNA_exprSet <- mRNA_exprSet[!index,]
row.names(mRNA_exprSet) <- mRNA_exprSet$gene_name
mRNA_exprSet$gene_name <- NULL


setwd('D:\\SCIwork\\F20ELFN1\\COAD')
save(mRNA_exprSet, file = "mRNA_exprSet.Rda")
第三步,提取肿瘤表达矩阵
##########################################################################################




###########################################################################################

rm(list=ls())


load( "mRNA_exprSet.Rda")



metadata <- data.frame(colnames(mRNA_exprSet ))

colnames(metadata) <- 'barcode'


for (i in 1:length(metadata[,1])) {
  num <- as.numeric(as.character(substring(metadata[i,'barcode'],14,15)))
  if (num == 1 ) {metadata[i,2] <- "T"}
  if (num != 1) {metadata[i,2] <- "N"}
}



names(metadata)[2] <- 'Barcode'

table(metadata$Barcode)

metadata <- subset(metadata, metadata$Barcode == 'T')


mRNA_exprSet <- mRNA_exprSet[,which(colnames(mRNA_exprSet) %in% metadata$barcode)]


setwd('D:\\SCIwork\\F20ELFN1\\COAD')

save(mRNA_exprSet, file = "mRNA_exprSet.Rda")

第四步,根据基因表达量筛选一些基因

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 一文解决大量基因的生存分析并作图

    这两篇纯生信文章都是对单个基因或者所有单个marker做生存分析,目的是找到其中能够影响患者生存的marker或者基因(包括miRNA,lncRNA,mRNA等...

    用户1359560
  • 重复一篇3分左右纯生信文章(第三部分)

    用户1359560
  • 单基因生信分析流程(1)一文解决TCGA数据下载整理问题

    在平常科研工作中,经常有师兄师姐师弟师妹问我:我现在有一个单基因,我该怎么开展生信研究?出现这个问题的原因是:(1)目前生信研究火热也逐渐受到认可(2)许多医学...

    用户1359560
  • 图解《移动终端支付可信环境技术规范》

    前面我们知道中国支付清算协会将TEE、SE作为一种产品类别,并将《移动终端支付可信环境技术规范》作为其技术规范。因为准确、全面的理解规范非常重要,本篇安智客将所...

    安智客
  • 为什么说Spark SQL远远超越了MPP SQLSpark SQL 成为了一种跨越领域的交互形态

    这里说的并不是性能,因为我没尝试对比过(下文会有简单的说明),而是尝试从某种更高一层次的的角度去看,为什么Spark SQL 是远远超越MPP SQL的。

    用户2936994
  • Typecho接入熊掌号附带AMP/MIP简单教程

    在Typecho接入熊掌号之前,你得先去熊掌号申请开通一下。具体填写的信息就是上面页面的那么多,请准备号正面手持身份证的照片。

    Erwin
  • WordPress设置页面的加载机制

    版权声明:本文为博主原创文章,遵循 CC 4.0 by-sa 版权协议,转载请附上原文出处链接和本声明。

    Jerry Wang
  • JVM解读-类加载机制

    Java虚拟机把描述类的数据从Class文件加载到内存,并对数据进行校验、转换解析和初始化,最终形成可以被虚拟机直接使用的Java类型,这就是虚拟机的加载机制。

    高广超
  • 你真的懂「类的加载机制」吗?

    高广超 :多年一线互联网研发与架构设计经验,擅长设计与落地高可用、高性能互联网架构。目前就职于美团网,负责核心业务研发工作。本文首发在 高广超的简书博客,欢迎点...

    用户1093975
  • 基于SSM框架搭建的项目,带你剖析MVC结构

    “ 这是小的Demo是我部署用来教大家学MVC小例子的,搭建完成SSM框架,非常简单,使用了Spring/Spring MVC/MyBatis框架,数据库使用了...

    赵腰静

扫码关注云+社区

领取腾讯云代金券