美文网首页
2018-09-06

2018-09-06

作者: 一路向前_莫问前程_前程似锦 | 来源:发表于2018-09-06 16:58 被阅读33次

    提起lncRNA

    options(stringsAsFactors = F)
    expr_df_nopoint=read.table("RNA_count_data.txt", sep = "\t",header = T)
    names(expr_df_nopoint)[1]="gene_id"

    首先我定义了一个非编码RNA的集合,这个每个人的标准不一样,但是我的原则是,多多益善,这样出来以后会有个问题,就是编码基因转录出非编码基因会无法从基因名称上区分,可以在运行时把geneid换成转录本id,必须要记在心里。

    详见:http://www.bioinfo-scrounger.com/archives/323

    lncRNA <- c("3prime_overlapping_ncRNA","antisense_RNA","bidirectional_promoter_lncRNA","lincRNA","macro_lncRNA","non_coding","processed_transcript","sense_intronic","sense_overlapping")

    读入GTF数据

    library(dplyr)
    library(rtracklayer)
    AnnoData <- rtracklayer::import('E:/GTF_R/Homo_sapiens.GRCh38.90.chr.gtf')
    AnnoData1 =data.frame(AnnoData )

    用 gene_biotype来筛选lncRNA

    LncRNA_exprSet <- AnnoData1 %>%
    filter(gene_biotype %in% lncRNA) %>% #注意这里是gene_biotype
    select(c(gene_name,gene_id,gene_biotype)) %>%
    distinct() %>% #删除多余行
    inner_join(expr_df_nopoint,by ="gene_id")%>%
    select(-gene_id,-gene_biotype)

    LncRNA_exprSet2<-LncRNA_exprSet%>%
    distinct(LncRNA_exprSetgene_name,.keep_all = T)%>% select(-`LncRNA_exprSetgene_name`)

    LncRNA_exprSet2`LncRNA_exprSetgene_name`

    Metadata 分组文件

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

    write.table(LncRNA_exprSet2, "LNC_RNA_count_data.txt", row.names = F, sep = "\t", quote = F)

    metadata <- data.frame(names(LncRNA_exprSet2)[-1])
    for (i in 1:length(metadata[,1]))
    {
    num <- as.numeric(substr(metadata[i,1],14,15))
    if (num %in% seq(1,9)) {metadata[i,2] <- "T"}
    if (num %in% seq(10,29)) {metadata[i,2] <- "N"}
    }

    加个名称

    names(metadata) <- c("TCGA_id","sample")

    将sample转化成factor

    metadatasample <- as.factor(metadatasample)
    table(metadatasample) class(metadatasample)

    我们可总结一下,有44例normal,502例肿瘤组织

    metadata %>% dplyr::group_by(sample) %>% summarise(n())

    ############# 说明所有的样本都是肿瘤样本

    制作countData

    Lnc_RNA_mycounts <- LncRNA_exprSet2

    row.names(Lnc_RNA_mycounts)<-Lnc_RNA_mycounts$gene_name
    Lnc_RNA_mycounts=Lnc_RNA_mycounts[-1]
    str(Lnc_RNA_mycounts)

    ############## log2(count+1) 归一化 #############
    tt=head(Lnc_RNA_mycounts)

    my_normalize<-function(x){
    x=log2(x+1)

    }

    tt2=apply(tt,2,my_normalize)

    Lnc_RNA_mycounts_norm<-apply(Lnc_RNA_mycounts,2,my_normalize)

    save(Lnc_RNA_mycounts_norm,file = "lncRNA_norm.Rda")

    相关文章

      网友评论

          本文标题:2018-09-06

          本文链接:https://www.haomeiwen.com/subject/pfgxgftx.html