RT
#加载数据 linba代表淋巴lnc <- read.csv("K0.csv")#保证ID的唯一性，删除重复table(duplicated(lnc$id))index<-duplicated(lnc$id)us_count<-lnc[!index,] #反向筛选#转化成矩阵，rowname转化成第一列A<-us_countrownames(A)<-A[,1] A<-A[,-1] #Count with fewer than 10 reads were removed.A[A < 10] <- 0#将输入数据取整A<-round(A,digits=0) A#准备，将数据转换为矩阵格式A<-as.matrix(A) ## 设置分组信息，建立环境(个样本，组处理)colnames(A)dim(A)condition <- factor(c(rep("control",179),rep("case",80)))coldata<-data.frame(row.names=colnames(A),condition) coldata #展示coldata内容condition #展示环境# 然后判断现在表达矩阵与样本信息是否一致，特别是针对很多样本时。all(rownames(coldata) == colnames(A))library("DESeq2") #使用library函数加载DEseq2包##构建dds矩阵 dds<-DESeqDataSetFromMatrix(A,coldata,design=~condition)head(dds) #去除表达量非常少的行，比如设置阈值为每行表达量为10。#基因至少在某一些文库的count超过10 ~ 15 才被认为是表达。keep <- rowSums(counts(dds)) >= 2000dds <- dds[keep,]dim(dds) #看过滤了#之前不设置时，默认control在case后面。 设置后，#control组(untreated)排在了前面#dds$condition <- factor(dds$condition, levels = c("control","case"))dds$condition <- relevel(dds$condition, ref = "control")levels(dds$condition)##进行差异分析#对原始的dds进行标准化dds<-DESeq(dds) #查看结果名称resultsNames(dds) #用results函数提取结果，并赋值给res变量res<-results(dds) #查看结果plotMA(res,ylim=c(-2,2)) mcols(res,use.names=TRUE)#简单绘制火山图plot(res$log2FoldChange,res$pvalue) #提取差异基因res <- res[order(res$padj),]resdata<-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)deseq_res<-data.frame(resdata)#提取上调差异表达基因up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 2)) #提取下调差异表达基因down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -2)) #输出上调基因write.csv(up_diff_result,"linbaUP.csv") #输出下调基因write.csv(down_diff_result,"linbaDOWN.csv") #输出总的差异write.csv(deseq_res,"linbaUpDown.csv")#简单验证下方向问题，看前面设置的对照组问题，#和上面导出的数据癌和对照的方向有没有错library(tidyverse)train <- up_diff_result[,-c(1:7)]train <- t(train)dim(train)train <- data.frame(train)train$group <-'case' dim(train)train$group[1:179] <- 'control'train$group <- factor(train$group)#通过平均值验证方向。#aggregate(x, by, FUN)aggregate(train[,1],list(train[,19]),mean)##火山图#提取差异的数据resdata <- mutate(resdata,select_change=if_else(abs(log2FoldChange) >= 1,'y','n'),select_pvalue=if_else( padj < 0.05,'y','n') )resdata <- mutate(resdata,select=paste(resdata[,'select_change'], resdata[,'select_pvalue'], sep = ''))table(resdata$select)library("ggplot2")##ggplot2 差异火山图resdata$select <- factor(resdata$select, levels = c('nn', 'ny', 'yn', 'yy'), labels = c('p >= 0.05, FC < 2', 'p < 0.05, FC < 2', 'p >= 0.05, FC >= 2', 'p < 0.05, FC >= 2'))#纵轴为显著性 p 值volcano_plot_pvalue <- ggplot(resdata, aes(log2FoldChange, -log(padj, 10))) +geom_point(aes(color = select), alpha = 0.6) +scale_color_manual(values = c('gray30', 'green4', 'blue2','red2')) +theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent')) +geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.5) + geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.5) +labs(x = 'log2 Fold Change', y = '-log10 p-value')volcano_plot_pvalue#保存为PDF#ggsave('volcano_plot_pvalue.pdf', volcano_plot_pvalue, width = 6, height = 5)
