前言

作为自己练习转录组的第一个项目，身为一个初学者，R使用的也不是特别熟练，估计会出现很多错误，而且一些处理方法我自己也会觉得很蠢，但是也没办法，只能一点点进步了。

分析步骤

表达矩阵的获取

我是利用生信人的TCGA下载小工具下载的LUSC数据，没有使用R包。最后获得了551例count数据，包括49例癌旁组织。

表达矩阵的预处理

library(stringr)
clinical = read.table('clinical.csv',header = T,sep = ',',row.names = 1)
names(clinical)
# [1] "A18_Sex"                             "A1_OS"                               "A22_histological_type"              
# [4] "A23_Tumor_tissue_site"               "A2_Event"                            "A3_T"                               
# [7] "A4_N"                                "A5_M"                                "A6_Stage"                           
#[10] "A8_New_Event"                        "A8_New_Event_Time"                   "A8_New_Event_Type"                  
#[13] "A9_Cancer_Status"                    "age_at_initial_pathologic_diagnosis" "anatomic_neoplasm_subdivision"
##这是我保留的我认为有用的临床信息。
rownames(clinical) = gsub('-','.',rownames(clinical)) 
##count的矩阵的分隔符是'.'，而临床信息的分隔符是'-'，先把它改成'.'
rt = read.table('counts.csv',header = T,row.names = 1,sep = ',')
rownames(rt) = str_split(rownames(rt),'[.]',simplify = T)[,1]
##去除基因名的版本号
clinical = clinical[clinical$A1_OS > 60,]
##筛选临床信息OS大于60天的病例，还有一些信息不完整的之前在excel中手动删除了，剩下456例样本
data = rt[,substr(names(rt),14,15) < 10]
normal = rt[,substr(names(rt),14,15) > 10]
##拆分癌症组织与癌旁组织
names(data) = gsub('.01$','',names(data))
##去掉.01
data = data[,!grepl('_Rep',names(data))]
names(data)
table(names(data)%in%rownames(clinical))
table(rownames(clinical)%in%names(data))
clinical = clinical[rownames(clinical)%in%names(data),]
data = data[,names(data)%in%rownames(clinical)]
##剩余453例样本
save(clinical,data,normal,rt,file = 'input.Rdata')

随后进行数据的预筛选与基因的注释

load('/path/input.Rdata')
library(edgeR)
library(stringr)
stopifnot(identical(rownames(data),rownames(normal)))
##为癌症样本结尾添加.01
names(data)  = gsub('(.*)','\\1.01',names(data))
##合并癌症癌旁组织
exprset = cbind(normal,data)
##分组
group_list=ifelse(as.numeric(substr(colnames(exprset),14,15)) < 10,'tumor','normal')
gene.names = rownames(exprset)
##添加基因注释信息
library(biomaRt)
mart = useMart('ensembl')
tmp = listDatasets(mart)
ensembl <- useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset="hsapiens_gene_ensembl",
                   host="aug2017.archive.ensembl.org") # Ensembl version 90.
ensemblGenes <- getBM(attributes=c('ensembl_gene_id',  'chromosome_name', 'gene_biotype', 
                                   'external_gene_name', 'entrezgene'), filters='ensembl_gene_id', 
                      values=gene.names, mart=ensembl) 
saveRDS(ensemblGenes, file="hg38.ensGene")
features <- ensemblGenes[match(gene.names, ensemblGenes$ensembl_gene_id),]
features$ensembl_gene_id <- gene.names
features$entrezgene <- gsub(" ", "", as.character(features$entrezgene))
row.names(features) <- gene.names
#features$Length <- gene.length
head(features)
##利用edgeR进行差异分析
y = DGEList(exprset,group = group_list)
y$samples
##cpm标准化
cpm = cpm(y)
lcpm = cpm(y,log = T)
##首先利用层次聚类查看离群样本
clust = t(lcpm)
clust = dist(clust)
pdf(file = 'hclust.pdf',width = 75)
plot(hc)
par(cex=0.01)
axis(1)
dev.off()
发现了这几个离群样本，在这里将他们删去
exprset = exprset[,!names(exprset) == 'TCGA.56.8623.01']
exprset = exprset[,!names(exprset) == 'TCGA.85.A4PA.01']
exprset = exprset[,!names(exprset) == 'TCGA.21.5782.01']
exprset = exprset[,!names(exprset) == 'TCGA.37.4129.01']
data = data[,!names(data) == 'TCGA.56.8623.01']
data = data[,!names(data) == 'TCGA.85.A4PA.01']
data = data[,!names(data) == 'TCGA.21.5782.01']
data = data[,!names(data) == 'TCGA.37.4129.01']
clinical = clinical[!rownames(clinical) =='TCGA.56.8623',]
clinical = clinical[!rownames(clinical) =='TCGA.85.A4PA',]
clinical = clinical[!rownames(clinical) =='TCGA.21.5782',]
clinical = clinical[!rownames(clinical) =='TCGA.37.4129',]
save(clinical,data,normal,exprset,file = 'input.Rdata')

#筛选lncRNA
ncRNA = c('3prime_overlapping_ncRNA',
'antisense_RNA',
'bidirectional_promoter_lncRNA',
'lincRNA',
'macro_lncRNA',
'non_coding',
'processed_transcript',
'sense_intronic',
'sense_overlapping')
LncRNA_exprSet<-features %>% dplyr::filter(features$gene_biotype %in% ncRNA)
lnc_expr = exprset[rownames(exprset)%in%LncRNA_exprSet$ensembl_gene_id,]

随后进行差异分析

#过滤基因
cpm = cpm(y)
lcpm = cpm(y,log = T)
keep <- rowMeans(lcpm>0)>0
y <- y[keep ,, keep.lib.sizes=FALSE]
dim(y)
#差异分析
y <- calcNormFactors(y, method = "TMM")
y$samples$norm.factors
#MDS图
cols = as.numeric(y$samples$group)
plotMDS(lcpm, col = cols)
title(main="Sample groups")

MDS

#构建design矩阵
design = model.matrix(~0+group_list)
colnames(design) <- gsub("group_list", "", colnames(design))
rownames(design) <- colnames(y)
design
#构建广义线性模型
par(mfrow=c(1,2))
y <- estimateDisp(y,design = design)
plotBCV(y)
y$common.dispersion
#只能说TCGA的样本量太大太杂，BCV已经是1.4了
fit <- glmQLFit(y, design, robust=TRUE)
plotQLDisp(fit)
con <- makeContrasts(tumor- normal, levels = design)
con
res <- glmQLFTest(fit, contrast=con)
summary(decideTestsDGE(res))

#火山图
plotMD(res)
DEmarkers <- topTags(res, n=Inf)$table
head(DEmarkers)
is.sig <- DEmarkers$FDR <= 0.05
plotSmear(res, de.tags=rownames(DEmarkers)[is.sig], cex=0.1)
#提取差异基因
diffSig = DEmarkers[(DEmarkers$FDR < 0.05 & (DEmarkers$logFC>2 | DEmarkers$logFC<(-2))),]#筛选有显著差异的#
write.table(diffSig, file="diffSig.txt",sep="\t",quote=F)#输出有显著差异表达的到diffSig这个文件#
diffUp = DEmarkers[(DEmarkers$FDR < 0.05 & (DEmarkers$logFC>2)),]#foldchange>0是上调，foldchange<0是下调#
write.table(diffUp, file="up.txt",sep="\t",quote=F)#39-42把上调和下调分别输入up和down两个文件#
diffDown = DEmarkers[(DEmarkers$FDR < 0.05 & (DEmarkers$logFC<(-2))),]
write.table(diffDown, file="down.txt",sep="\t",quote=F)