读取基因对应类型的数据
rm(list=ls())
## for human gencode.v25.annotation.gtf
# cat gencode.v25.annotation.gtf|perl -alne '{next unless $F[2] eq "gene";print}'|grep -w HAVANA |\
# cut -f 1,4,5,9| cut -d";" -f 1,2,4|sed 's/gene_id//g'|sed 's/gene_type//g'|sed 's/gene_name//g'|\
# sed 's/;//g'| sed 's/\"//g'|perl -alne '{/(ENSG\d+)/;print "$1\t$_"}' >human.gene.positions
if(F){
a=read.table('data/human.gene.positions')[,c(2:4,1,6,7)]
colnames(a)=c('chr','start','end','ensembl','type','symbol')
length(unique(a$symbol))
length(unique(a$ensembl))
head(a)
human_geneInfo_genecode_v25=a
}#生成基因及类型对应表，见下图
QQ截图20190714223847.jpg
load("D:/R/tcga_example-master/tcga_example-master/human_geneInfo_genecode_v25.rda")
#load('../human_geneInfo_genecode_v25.rda')
b=human_geneInfo_genecode_v25
head(b)
tail(sort(table(b$type)))#按照类型排序
QQ截图20190714224348.jpg
#lincRNA_genes=b[b$type %in%'lincRNA','ensembl']#选lincRNA对应的基因，基因名取的ensembl
lincRNA_genes=b[b$type %in% c('lincRNA','antisense'),'ensembl']
#现在定义的lincRNA就包括antisense和lincRNA
#%in%挑选 b$type里面属于 c('lincRNA','antisense')的出来
读取表达矩阵
a=read.table('Rdata/TCGA-UVM.htseq_counts.tsv.gz',header = T,sep = '\t')
dim(a)
a[1:4,1:4]
把表达矩阵的中对应的lincRNA的基因挑选出来
library(stringr)#通过字符串函数对字符串进行处理
a$Ensembl_ID=str_split(a$Ensembl_ID,'[.]',simplify = T)[,1]
#str_split处理字符串，用点进行分割成两部分后取前面的部分
lincRNA=a[a$Ensembl_ID %in% lincRNA_genes,]
#再用 %in% 进行挑选
#后续可以用lincRNA做下游分析
