RNA-seq的數(shù)據(jù)分析是比較簡(jiǎn)單基礎(chǔ)的分析，大概流程就是處理下機(jī)的fastq數(shù)據(jù)（trimmomatic）忱屑，比對(duì)到人類基因組（hisat2）然后統(tǒng)計(jì)每個(gè)基因上出現(xiàn)的counts數(shù)（featureCounts），接下來(lái)在R里進(jìn)行差異表達(dá)分析（DEseq2）找出差異表達(dá)基因再進(jìn)行一些富集分析（clusterprofiler）函荣。
因?yàn)榍皫滋靹偤锰幚砹艘慌?0個(gè)樣本的RNA-seq數(shù)據(jù)绢记，我把每一步都記錄下來(lái)令杈。
首先是在linux下進(jìn)行處理，得到每個(gè)樣本的counts數(shù)文件匀借。

#!/bin/sh

fq_path=/work/data/lch_analysis
out_path="/work/analysis/lch_analysis"
trim_path="/opt/Trimmomatic-0.35"
hisat2_path="/opt/"
ht2_genome="/work/database/human/HISAT2/UCSC_hg19_Human"
gtf="/work/database/human/GTF/gencode.v27lift37.annotation.gtf"

cd $fq_path || exit 1
for file in $(ls | grep _1.fq.gz)
        do
        pre_name=${file%_good_1.fq.gz}
        mkdir -p $out_path/$pre_name/log
        fq1="${pre_name}_good_1.fq.gz"
        fq2="${pre_name}_good_2.fq.gz"

        # Step 0
        # Trim reads using Trimmomatic
        java -Xmx4g -jar \
        $trim_path/trimmomatic-0.35.jar \
        PE \
        -threads 4 \
        -phred33 \
        $fq1 \
        $fq2 \
        -baseout $out_path/$pre_name/${pre_name}_trimmed.fq.gz \
        ILLUMINACLIP:$trim_path/adapters/TruSeq2-PE.fa:2:30:10:6:true \
        SLIDINGWINDOW:4:15 \
        MINLEN:51 \
                2> $out_path/$pre_name/log/${pre_name}_Trimmomatic.log

        # Step 1
        # Map reads to hg19 reference
        hisat2 \
        -p 4 -q \
        -x $ht2_genome/genome \
        --fr \
        --rg-id ${pre_name} \
        --rg SM:${pre_name} --rg LB:${pre_name} --rg PL:ILLUMINA \
        -1 $out_path/$pre_name/${pre_name}_trimmed_1P.fq.gz \
        -2 $out_path/$pre_name/${pre_name}_trimmed_2P.fq.gz | samtools sort -@ 4 - \
        -T $out_path/$pre_name/log/rna_temp \
        -l 1 \
        -o $out_path/$pre_name/${pre_name}_sorted.bam \
        2> $out_path/$pre_name/log/${pre_name}_samtools_sort.log

        # Step 2
        # Count number of reads on genes
        featureCounts \
                -T 4 \
                -p \
                -t exon \
                -g gene_id \
                -a $gtf \
                -o $out_path/$pre_name/${pre_name}_featureCounts.txt \
                $out_path/$pre_name/${pre_name}_sorted.bam \
                2> $out_path/$pre_name/log/${pre_name}_featureCounts.log

        done

上述步驟挺簡(jiǎn)單的颜阐，就不贅述了。接下來(lái)是在R中的處理吓肋，因?yàn)楣灿?0個(gè)文件凳怨，所以要批量讀入處理，這里就要用到lapply和sapply函數(shù)了蓬坡。
lapply的使用格式為：
lapply(X, FUN, ...)
lapply的返回值是和一個(gè)和X有相同的長(zhǎng)度的list對(duì)象猿棉，這個(gè)list對(duì)象中的每個(gè)元素是將函數(shù)FUN應(yīng)用到X的每一個(gè)元素磅叛。其中X為L(zhǎng)ist對(duì)象（該list的每個(gè)元素都是一個(gè)向量）屑咳，其他類型的對(duì)象會(huì)被R通過(guò)函數(shù)as.list()自動(dòng)轉(zhuǎn)換為list類型。
函數(shù)sapply是函數(shù)lapply的一個(gè)特殊情形弊琴，對(duì)一些參數(shù)的值進(jìn)行了一些限定兆龙，其使用格式為：
sapply(X, FUN,..., simplify = TRUE, USE.NAMES = TRUE)
sapply(x, simplify = FALSE, USE.NAMES = FALSE) 和lapply()的返回值是相同的。如果參數(shù)simplify=TRUE敲董，則函數(shù)sapply的返回值不是一個(gè)list紫皇，而是一個(gè)矩陣；若simplify=FALSE腋寨，則函數(shù)sapply的返回值仍然是一個(gè)list聪铺。

source("http://bioconductor.org/biocLite.R")
biocLite("DESeq2")
library(DESeq2)
library(org.Hs.eg.db)
library(pheatmap)
setwd('/work/work/rna-seq-lhc/deseq2_results')
##讀取同一目錄下的所有文件
path <- "/work/work/rna-seq-lhc/featurecounts_results_cut" ##文件目錄
fileNames <- dir(path)  ##獲取該路徑下的文件名
filePath <- sapply(fileNames, function(x){ 
  paste(path,x,sep='/')})   ##生成讀取文件路徑
data <- lapply(filePath, function(x){
  read.table(x,sep = '\t', header=F,stringsAsFactors = F)})  ##讀取數(shù)據(jù)，結(jié)果為list

data2 <- lapply(data, function(x){data.frame(ID=gsub('\\..*', '', x$V1), Count=x$V2)})##去掉ensemble版本號(hào)
data3 <- do.call(cbind, data2)##把60個(gè)樣本合并
data4 <- data3
data4$ID <- data4$`Human_J429-ZX01-L01_featureCounts.txt.ID`
data4 <- data4[,c('ID', grep('Count$', colnames(data4), value = T))]
data5=aggregate(data4[,-1],by=list(data4$ID),sum)##重復(fù)的基因相加
rownames(data5)=data5$Group.1
data5=data5[,-1]
colnames(data5)=unlist(strsplit(colnames(data5),'_f'))[seq(1,ncol(data5)*2,by=2)]
data5[1:6,1:6]

table=read.csv('/work/work/rna-seq-lhc/分組信息.csv')
table1=table[table$condition=='before'|table$condition=='after',]##選出治療前和治療后的病人
count1=data5[,table1$name]
countmatrix1<-as.matrix(count1)
dds1 <- DESeqDataSetFromMatrix(countmatrix1, colData=table1, design= ~ condition)
dds1 <- dds1[ rowSums(counts(dds1)) > 1, ]
dds1 <- DESeq(dds1)##標(biāo)準(zhǔn)化
res1 <- results(dds1)
res1$genesymbol <- mapIds(org.Hs.eg.db,keys = rownames(res1),column = "SYMBOL",keytype ="ENSEMBL",multiVals = 'first')##加一列基因symbol
res1 <- res1[order(res1$padj),]
res1 <- res1[,c(7,1,2,3,4,5,6)]
res1 <- merge (as.data.frame(res1),as.data.frame(counts(dds1,normalize=TRUE)),by="row.names",sort=FALSE)
write.table(res1[1:1000,],"result_before_after_1000.csv", sep = ",", row.names = T)
differgenes1<-subset(res1,padj<0.05&(log2FoldChange>1|log2FoldChange< -1))

log2.norm.counts1 <- as.data.frame(log2(counts(dds1,normalize=T)+1))[differgenes1$Row.names,]
log2.norm.counts1[log2.norm.counts1>3]=3
log2.norm.counts1 <- t(scale(t(log2.norm.counts1)))##scale
annotation_col1 <- data.frame(condition=table1$condition)
rownames(annotation_col1) <- table1$name
pheatmap(log2.norm.counts1, cluster_rows=TRUE, show_rownames=FALSE,
         cluster_cols=T, annotation_col = annotation_col1)

輸出的res1就是差異表達(dá)前一千的基因萄窜，后面是DEseq2標(biāo)準(zhǔn)化后的counts數(shù)铃剔。如果有需要，還可以做一些tsne圖查刻，火山圖键兜，富集圖。我并沒(méi)有做穗泵，在這里把之前相關(guān)的腳本貼一下普气。
火山圖

#valcano plot
library(ggplot2)
library(openxlsx)
library(ggrepel)
library(dplyr)
data=read.xlsx('~/桌面/初步分析2-egfl7.xlsx',sheet = 1)
#data=read.xlsx('~/桌面/初步分析-mir126.xlsx',sheet = 1)

data$threshold <- as.factor(ifelse(data$pvalue < 0.05 & abs(data$log2FC) >= 1,ifelse(data$log2FC > 1 ,'Up','Down'),'Not'))
p=ggplot(data=data,aes(x=log2FC, y =-log10(pvalue),colour=threshold,fill=threshold)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_point(alpha=0.8, size=1.2)+
  xlim(c(-4, 4)) +
  theme_bw(base_size = 12, base_family = "Times") +
  geom_vline(xintercept=c(-1,1),lty=4,col="grey",lwd=0.6)+
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)+
  theme(legend.position="right",
        panel.grid=element_blank(),
        legend.title = element_blank(),
        legend.text= element_text(face="bold", color="black",family = "Times", size=8),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(face="bold", color="black", size=12),
        axis.text.y = element_text(face="bold",  color="black", size=12),
        axis.title.x = element_text(face="bold", color="black", size=12),
        axis.title.y = element_text(face="bold",color="black", size=12))+
  labs(x="log2 (fold change)",y="-log10 (p-value)",title="Volcano picture of DEG")

p+geom_text_repel(data=filter(data, pvalue< 4.32E-09), aes(label=genesymbol),show_guide=F)##把p值小于4.32E-09的基因標(biāo)注出來(lái)

tsne圖

##tsne
tsne_matrix <- t(degs_counts)
tsne_result <- Rtsne(tsne_matrix,perplexity = 3,pca = F,theta=0.5)
tsne_plot <- data.frame(Cluster.1 = tsne_result$Y[,1], Cluster.2 = tsne_result$Y[,2], 
                        Type = factor(annotation$Type,levels = c("health","before")))
ggplot(tsne_plot) + 
  geom_point(aes(x=Cluster.1, y=Cluster.2, color=Type)) +
  theme_bw(base_size = 12, base_family = "") +
  theme(legend.justification=c(0,0),legend.position = c(0.75,0.75),
        legend.title = element_blank(),
        legend.text = element_text(size = 10)) + 
  ggtitle("t-SNE Clustering (top 100 DEGs)") + 
  theme(plot.title = element_text(hjust = 0.5))
dev.off()

GO、KEGG富集

source("https://bioconductor.org/biocLite.R")
biocLite("org.Hs.eg.db")
setwd("/work/R語(yǔ)言/")
library(org.Hs.eg.db)
#install clusterProfiler
source("https://bioconductor.org/biocLite.R")
biocLite("org.Mm.eg.db")
library(org.Mm.eg.db)
biocLite("clusterProfiler")
library(clusterProfiler)
#GO
ego<-enrichGO(OrgDb="org.Mm.eg.db", 
             #gene = row.names(differgenes),
             gene = rownames(results[grep("1",results$GeneCluster),]),
             pvalueCutoff = 0.01,
             keytype = "ENSEMBL",
             readable=TRUE)
write.csv(as.data.frame(ego),"G-enrich.csv",row.names =F)

#KEGG
a=read.csv("/work/new_output/filter_immune/at2_subtype/at2.csv",header = T,sep = ",",stringsAsFactors = F)
x<-select(org.Mm.eg.db,
          keys = a[,z], 
          column = "ENTREZID", 
          keytype = "ENSEMBL"
          )
kegg<-x10[,2]
ekk <- enrichKEGG(gene=kegg,
                  keyType = "kegg",
                  organism = 'mmu',
                  pvalueCutoff = 0.05,
                  pAdjustMethod = "BH", 
                  qvalueCutoff = 0.1)
DOSE::dotplot(ekk, font.size=10)
  write.csv(as.data.frame(ekk),y,row.names =F)