NPC-103611中的id转换问题

关于NPC(鼻咽癌)一个数据集GSE103611的前期探索，暂且记录一下。一开始发现这个gene_assignment这个列很复杂。

image-20200610130030671

总结就是NM_是蛋白质编码RNA，NR_是非编码RNA,也就是如下

image-20191114113950770

image-20191114114150213

然后去那个npc那个数据集去看看大小

mrna <- ids[grep('NM',ids$gene_assignment),]
nr_rna <- ids[grep('NR',ids$gene_assignment),]

#从nr_rna里找还带有NM的,发现有NM
a <- grep('NR',nr_rna$gene_assignment)
b <- grep('NM',nr_rna$gene_assignment)
intersect(a,b)

#从mrna里找还带有NR的,发现有NR
c <- grep('NM',mrna$gene_assignment)
d <- grep('NR',mrna$gene_assignment)

得到的mrna和nr_rna就是分别对应的提取出的按照NM和NR,得到大小如下

image-20191114114914699

image-20191114114322044

上面这张图显示着总共有14265行是既含有NM也含有NR的，那没关系我去掉就好了，即使去掉的话NR也还有16853-14265个，还有1万多个,而对于蛋白编码的NM来说，那去掉1万多个还有19万多个呢

==1.我的第一个问题是，能这么去掉嘛？那些既含有NR又含有NM的到底是蛋白编码基因还是非编码基因呢？如果不想弄懂，就都去掉也行，但不知道要不要这么做，会不会丢失信息==

==2.然后第二个问题是，老大，是在我想问你时想到的，就是本疑问蛋白编码有近20万个，做WGCNA电脑会奔溃，因为在wgcna的tutois上到说好像50000个探针算很高了==然后我想到上午我点那个==nr_rna====发现那个其实是很多NM对应着一个基因的，那么把重复的去掉就好了，这样就不会是近20万个探针了，然后接下就做个事==

image-20191114113705536

==这个我也想不通下面这个是两个按照NR筛选出的nr_rna的，但是出现了NM号，后面还跟着NR号而且连续挨着的，一上一下两个，他们的细节都一样，只有ENST的顺序不一样，为什么还列出来了呢==，那我就把这样的都去掉吧老大

image-20200610131322927

image-20200610131342072

==老大说，用genecode注释==

现在该如何做起？思路如下：

1.现在切除来的那列是都是基因名，但是基因名都是重复的，所以我应该先不着急区分出是编码还是非编码

2.应该是把重复的基因名，merge到表达矩阵上，然后根据基因表达量，只留表达量最大的探针，或者是取中位值

3.等基因名全部都是uniqe的时候，通过genecode来进行注释，分成ncRNA和mRNA.

image-20191114180847565

linux里输入如下,保存为gtf.txt的文件

zless -S gencode.v32.annotation.gtf.gz|grep -w 'gene'|cut -f 9|cut -d ';' -f 1-3|tr ';' ' '|sed 's/"//g'|awk '{print $4,$6 }'|less -S > gtf.txt

false的原因是有空格，gsub去掉

image-20191115152631523

转换行名时开始有报错

rownames(ex) <- gsub(' ','',rownames(ex))

image-20200610131410298

image-20191115162037423

这个地方很有问题，因为我的rownames都是前面根据median值给去重复的了，可是现在竟说有重复的,那我就新建一列。

rm(list = ls())
options(stringsAsFactors = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
library(GEOquery)
library(WGCNA)
#exprset_gse103611<-read.csv('GSE103611_series_matrix.txt',comment.char = '!',sep = '\t',header = T)
f<-'GSE103611.Rdata'
if(!file.exists(f)){
  gset<-getGEO('GSE103611',destdir='.',
               AnnotGPL=F,
               getGPL=F)
  save(gset,file=f)
}
load('GSE103611.Rdata')
ex<- exprs(gset[[1]])
dim(ex)
ex[1:4,1:4]
pd <- pData(gset[[1]])
pd <- pd[,apply(pd,2,function(x){
  length(unique(x))> 1
})]
colnames(pd)
datTraits<-data.frame(row.names=rownames(pd),metastasis_status=pd[,7])
dim(datTraits)

ids<-read.csv('GPL19251-14.txt',sep = '\t',comment.char = '#',header = T)
ex <- as.data.frame(ex)
save(ex,ids,file = 'step1_raw.Rdata')
load('step1_raw.Rdata')

ids[1:4,1:4]
head(ids)
ids <- ids[ids$ID %in% rownames(ex),] #此时用的%in% 并不会改变ids的顺序，ids是一个一个在ex里面实验，在就列上去，不在就扔掉，继续下一个
head(ids)
dim(ids)
# mrna <- ids[grep('NM',ids$gene_assignment),]
# nr_rna <- ids[grep('NR',ids$gene_assignment),]
# 
# #从nr_rna里找还带有NM的,发现有NM
# a <- grep('NR',nr_rna$gene_assignment)
# b <- grep('NM',nr_rna$gene_assignment)
# intersect(a,b)
# 
# #从mrna里找还带有NR的,发现有NR
# c <- grep('NM',mrna$gene_assignment)
# d <- grep('NR',mrna$gene_assignment)
# 
# ##上面说明那么就是筛选后的mrna和nr_rna里，是有那些既含有NR也含有NM的基因的。那么是什么基因呢

ids$gene <- str_split(ids$gene_assignment,'//',simplify = T)[,2]
id2symbol <- ids[,c(1,14)]
dim(id2symbol)
colnames(id2symbol)=c('probe_id','symbol')
id2symbol <- id2symbol[id2symbol$symbol !='',]
dim(id2symbol)

# id2symbol <- id2symbol[id2symbol$probe_id %in% rownames(ex),]
# dim(id2symbol)
ex[1:4,1:4]
ex <- as.data.frame(ex)
ex <- ex[match(id2symbol$probe_id,rownames(ex)),]
#ex <- ex[id2symbol$probe_id,]  ####这步因为探针名也纯是数字
head(rownames(ex))
dim(ex)
dim(id2symbol)
ex[1:4,1:4]



id2symbol$median=apply(ex,1,median)
head(id2symbol)
id2symbol=id2symbol[order(id2symbol$symbol,id2symbol$median,decreasing = T),]
id2symbol=id2symbol[!duplicated(id2symbol$symbol),]#
dim(id2symbol)
id2symbol$probe_id <- as.character(id2symbol$probe_id)

dim(ex)  #到现在位置，ex的探针顺序应该和id2symbol的probe_id这一列的顺是完全一致的，因为是按照行名取的子集
ex <- ex[match(id2symbol$probe_id,rownames(ex)),]
table(id2symbol$probe_id %in% rownames(ex))


save(ex,id2symbol,file = 'step1-1.Rdata')


gtf <- read.table('gtf.txt',sep = ' ')
head(gtf)
colnames(gtf) <- c('gene_type','symbol')
identical(id2symbol$probe_id,rownames(ex))
rownames(ex) <- id2symbol$symbol

rownames(ex) %in% gtf$symbol
table(rownames(ex) %in% gtf$symbol)
head(rownames(ex))
head(gtf$symbol)

ex$genename <- gsub(' ','',rownames(ex))
table(ex$genename %in% gtf$symbol)

x <- intersect(ex$genename,gtf$symbol)
exfil <- ex[match(x,ex$genename),]
gtf_fil <- gtf[match(x,gtf$symbol),]

identical(exfil$genename,gtf_fil$symbol)
table(gtf_fil$gene_type)

lnc_gtf <- gtf_fil[gtf_fil$gene_type=='lncRNA',]
mi_gtf <- gtf_fil[gtf_fil$gene_type=='miRNA',]
mrna_gtf <- gtf_fil[gtf_fil$gene_type=='protein_coding',]

lnc_ex <- exfil[match(lnc_gtf$symbol,exfil$genename),]
mi_ex <- exfil[match(mi_gtf$symbol,exfil$genename),]
mrna_ex <- exfil[match(mrna_gtf$symbol,exfil$genename),]
rownames(lnc_ex) <- lnc_ex$genename
rownames(mi_ex) <- mi_ex$genename
rownames(mrna_ex ) <- mrna_ex$genename
lnc_ex <- lnc_ex[,-49]
mi_ex <- mi_ex[,-49]
mrna_ex <- mrna_ex[,-49]



save(lnc_gtf,mi_gtf,mrna_gtf,lnc_ex,mi_ex,mrna_ex,file = 'step1-2.Rdata')

image-20200610131430539

这个是比较简单的事情，就是把’yes‘改成’NPC‘,'no'改成‘normal’。

image-20200610131441557

报错如下：==不知道怎么改了，就是说不是yes的有很多个，就是no的个数嘛，用一个‘normal’去赋值给给那么多的no就报错了，但是这个地方不知道怎么改==

image-20191115235134303

后面的分析现在看来都不对，首先下面这个层次聚类图就不对，17到40号应该都是转移组，二41到64应该是非转移组，现在二者完全分不开，混为一锅粥

image-20191116094901080

看一下层次聚类的概念，参考层次聚类算法的原理及实现Hierarchical Clustering:https://blog.csdn.net/zhangyonggang886/article/details/53510767

image-20200610131523221

上面的两张图和我内心的期待是一样的,一样的道理咯，就是表达量相当于距离。

==所以，我得到的那种聚类树的图应该是不正确的==

现在思考有没有可能是注释文件不对，我用的是gencode版本的，但是看下图，我发现这个红框圈住的基因是来自ensembl的，我想的是有没有==可能就是genecode的文件注释不全==，再下载一遍ensembl的gtf文件

image-20191116104940735

image-20200610131537393

没有重新下载，现在做一步法的时候，明明net$colors是有好多个的颜色的模块的，如下图

image-20191116160410611

但是下面的却图应该是给自动合并了，我不知道如何调参数，代码也是之前的代码

image-20200610131553274

找到了原因

  sizeGrWindow(12, 9)
  mergedColors = labels2colors(net$colors)
  pdf('step3-dynamicColors_mergedColors.pdf',width = 25,height =20 )
  # Plot the dendrogram and the module colors underneath
  plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                      "Module colors",
                      dendroLabels = FALSE, hang = 0.03,
                      addGuide = TRUE, guideHang = 0.05)

是因为这个plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],选择的是第一个基因分类，而net$dendrograms一共有10个list，如果将==1==改为2、3出图不一样，但是我不会将图片合并。之前之所以是==1==，是代码中我增加了maxBlockSize = 2000,那么就把1万7千多的探针给分成了每组不超过2000的基因块，这个macBlockSize默认是5000，而前面老大取了就是根据mad值取了5000个基因，这样取出来==[[1]]==是正好的。

image-20191116161453697