转录组全新学习之总篇

Y叔的RNA-seq workflow 强力推荐生信技能树讲解 1. 数据已经存在服务器，基因组注释文件、基因组文件都已经存在。服务器环境软件都已经安装完成。但是没有root权限。数据获取预处理

find ./510 -name '*zaoqian*'|xargs -i mv {} ./100    #批量查找510目录里所有包含zaoqian的文件，移动到100目录里。

参考 1 2

实验流程和使用软件
质控 trimmomatic trimmomatic功能讲解 精华篇 trimmomatic参数讲解
比对 HISAT2 精华HISAT2 HISAT2官网
装配 String tie 精华stringtie
定量分析 Deseq2
转录本分析 Cufflink
从头组装 trinity
开始第一步——质控。 fastqc可以检测测序质量。

fastqc *.fastq.gz -t 4           #开启4个进程。

fastqc测序质量评估

trimmomatic去除接头和低质量测序数据（-threads 12 开启了12个进程）

java -jar /home/guo/tool/Trimmomatic-0.38/trimmomatic-0.38.jar PE -threads 12 -phred33 CR-zaoqian-1_L8_I313.R1.clean.fastq.gz CR-zaoqian-1_L8_I313.R2.clean.fastq.gz CR-paired-1-R1.fastq.gz CR-unpaired-1-R1.fastq.gz CR-paired-1-R2.fastq.gz CR-unpaired-1-R2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 &

最后的&是在后台执行。 4. 比对HISAT2 - 先使用hisat2的自带的py脚本处理基因组注释文件和snp文件

#运行目录是 /disks/backup/chaim/maize/
/home/chaim/disk/soft/hisat2/extract_exons.py Zea_mays.B73_RefGen_v4.42.gtf > genome.exon
/home/chaim/disk/soft/hisat2/extract_splice_sites.py Zea_mays.B73_RefGen_v4.42.gtf > genome.ss
/home/chaim/disk/soft/hisat2/hisat2_extract_snps_haplotypes_VCF.py zea_mays.vcf> genome.snp

建立索引文件（此步骤非常耗时，至少1-2h,多则1-2天）

#-p 8 8个进程，基因组文件19v4.dna.fa  /disks/···/maize437 指定输出的文件的名字为maize437，存放位置就是前面这一串路径。
#hisat2-build -p 8 /disks/backup/chaim/maize/Zea_mays.B73_RefGen_v4.42.fa --snp genome.snp --ss genome.ss --exon genome.exon /disks/backup/chaim/maize/genome_tran &~~
#鉴于添加snp之后，服务器实在吃不消。所以就不添加snp位点了。内存在200G以上的可以试一试加snp的建立索引。
 hisat2-build -p 8 /disks/backup/chaim/maize/Zea_mays.B73_RefGen_v4.42.fa --ss genome.ss --exon genome.exon /disks/backup/chaim/maize/genome_tran &

开始比对


hisat2 -x /disks/backup/chaim/maize/genome_tran -p 8 -1 /disks/backup/chaim/cms/zaoqian/Cr-paired-1-R1.fastq.gz -2 /disks/backup/chaim/cms/zaoqian/Cr-paired-1-R2.fastq.gz -S Cr-1.map.sam --dta-cufflinks --novel-splicesite-outfile Cr-1.nsplice &

参数讲解 -x /disks/backup/chaim/maize/genome_tran #基因组索引文件(后缀是.fai,此处只写文件名，不要.fai后缀) -p 8 #使用8个进程 -1 paired1 #质控后左链的数据cms -2 paired2 #质控后右链的数据 -S Cr-1.map.sam #输出文件 -dta-cufflinks--novel-splicesite-outfile Cr-1.nsplice #后期使用cufflinks做准备。 4.5 samtools将sam转bam同时排序。

samtools sort -@ 8 -o ${bef[$i]}".map.bam" ${bef[$i]}".map.sam" 2>${bef[$i]}"samtool_out" 


#构建索引文件，方便使用IGV查看拼接结果
samtools index -@ 8 {bef[$i]}".map.bam" {bef[$i]}".map.bai" &

stringtie官网文档（自带梯子）

for (( i=1;i<4;i++ ));
do
stringtie CR.map.bam -G /disks/backup/chaim/maize/19v4.42.gff3 -p 8 -o CR-${i}.gtf &
stringtie Nr.map.bam -G /disks/backup/chaim/maize/19v4.42.gff3 -p 8 -o Nr-${i}.gtf &
stringtie Cr.map.bam -G /disks/backup/chaim/maize/19v4.42.gff3 -p 8 -o Cr-${i}.gtf &
done


stringtie --merge -G /disks/backup/chaim/maize/19v4.42.gff3 -o merged.gtf CR-1.gtf CR-2.gtf CR-3.gtf Cr-1.gtf Cr-2.gtf Cr-3.gtf Nr-1.gtf Nr-2.gtf Nr-3.gtf &

for (( i = 0; i < 4; i++ )); 
do
    stringtie Cr-${i}.map.bam -G merged.gtf -p 8 -b ${A}${i}"_out" -e -o Cr-${i}ss.gtf &
    stringtie Nr-${i}.map.bam -G merged.gtf -p 8 -b ${B}${i}"_out" -e -o Nr-${i}ss.gtf &
    stringtie CR-${i}.map.bam -G merged.gtf -p 8 -b ${C}${i}"_out" -e -o CR-${i}ss.gtf &
done

参数讲解： -G 参考基因组注释文件（我使用 gtf格式的注释，到第三次会出错，所以使用的是gff3的注释） -p 8个进程 -o 输出注释文件名 -b 输出其他文件的路径（会同时生成多个文件，一定要输出在不同的路径中，不然后面的输出会覆盖前面的输出。） 5.1 stringtie的其他功能——gffcompare和gffread 生信技能树参考官方地址

6.格式转换，数据整合。使用的是python的脚本。下载地址

python2.7 /disks/backup/chaim/soft/prepDE.py -i gtf2
#注意此处的prepDE.py的python版本是python2.7，一定要指定版本号。否则会报错。

gtf2文件内容是对应的比对后文件名和文件的位置。

Cr1-st ./Cr1-st.gtf
Cr2-st ./Cr2-st.gtf
Cr3-st ./Cr3-st.gtf
CR1-st ./CR1-st.gtf
CR2-st ./CR2-st.gtf
CR3-st ./CR3-st.gtf
Nr3-st ./Nr3-st.gtf
Nr1-st ./Nr1-st.gtf
Nr2-st ./Nr2-st.gtf

输出文件为转录本输出矩阵 transcript_count_matrix.csv 基因输出矩阵 gene_count_matrix.csv

perpDE.py参数

-i 输入文件 -g gene输出矩阵位置 -t 转录本输出矩阵位置

Deseq2定量分析参考地址 ~~安装在本地Ubuntu环境~~ ~~使用conda安装~~ ~~>conda install biocpnductor-deseq2~~

安装到服务器端 source("[http://bioconductor.org/biocLite.R](http://bioconductor.org/biocLite.R)") biocLite("DESeq2", dependencies=T)

使用DESeq2

library("DESeq2")
countData <- as.matrix(read.csv("gene_count_matrix.csv",row.names = "gene_id"))
head(countData,10)
#               CR1.st CR2.st CR3.st Cr1.st Cr2.st Cr3.st Nr1.st Nr2.st Nr3.st
#MSTRG.28369         0      0      0      0      0     68      0      0      0
#Zm00001d015574      8     12      0      0      4      0      8      0      0
#Zm00001d015575      0      0      5      7      0      0      0      9      0
#Zm00001d007250      0      0      0      0      0      0      0      0      0
#Zm00001d007257      0      0     12      0      0      0      0      0      0
#Zm00001d015572      0      0     17      0      0      0      0      9     12
#Zm00001d015573      0      0      0      0      0      0      0      0      0
#Zm00001d015578      0      0     13      0      0     58      0      0      0
#Zm00001d015579      0      4     23     39      0     10      0      0      5

condition <- factor(c(rep("CR",3),rep("Cr",3),rep("Nr",3)),levels=c("CR","Cr","Nr"))
###此处的levels向量，CR相当于control，其他的是以CR为基准比较上调或下调。deseq2默认的levels是c("未处理","处理")一定要注意这个顺序。不然可能会产生相反的上升和下降值。
colData <- data.frame(row.names = colnames(countData),condition)
head(colData,20)
#   condition
#CR1.st        CR
#CR2.st       CR
#CR3.st       CR
#Cr1.st        Cr
#Cr2.st        Cr
#Cr3.st        Cr
#Nr1.st        Nr
#Nr2.st        Nr
#Nr3.st        Nr

##上面构建colData和countData,应根据自己实际样本数略作调整。但应始终保持countData的第一行和colData的第一列完全一致，不包括countData的第一个空格。

dds <- DESeqDataSetFromMatrix(countData = countData,colData = colData, design = ~ condition)

dds <- DESeq(dds)
dds
#class: DESeqDataSet
#dim: 56862 9
#metadata(1): version
#assays(4): counts mu H cooks
#rownames(56862): MSTRG.28369 Zm00001d015574 ... #Zm00001d050294
 # Zm00001d005327
#rowData names(26): baseMean baseVar ... deviance maxCooks
#colnames(9): CR1.st CR2.st ... Nr2.st Nr3.st
#colData names(2): condition sizeFactor
res = results(dds)
res = res[order(res$pvalue),]
head(res)
#log2 fold change (MLE): condition Nr vs CR
#Wald test p-value: condition Nr vs CR
#DataFrame with 6 rows and 6 columns
  #                  baseMean   log2FoldChange             lfcSE
   #               <numeric>        <numeric>         <numeric>
#MSTRG.14834 813.537803375226 5.92709649671671 0.505759997283559
#MSTRG.39476 1429.58827897959 9.72303418884463 0.897070409158461
#MSTRG.14858 5517.88559178615 7.13535760976146 0.676567107148184
#MSTRG.14860 634.902462467296 8.24781803490373 0.955826315936781
#MSTRG.14664 1638.94705588096 7.01525377811048 0.823766922503449
#MSTRG.15005 1570.64339845137 9.18561797686647  1.12651942455926
  #                      stat               pvalue                 padj
  #                  <numeric>            <numeric>            <numeric>
#MSTRG.14834 11.7191880112132 1.01641526263862e-31 3.71276167136634e-27
#MSTRG.39476 10.8386522279403 2.25772589986083e-27 4.12351058350581e-23
#MSTRG.14858 10.5464151809536 5.27723078445638e-26 6.42555620315409e-22
#MSTRG.14860 8.62899241984173 6.18947935284045e-18  5.6522325450139e-14
#MSTRG.14664 8.51606636109028 1.65064269019706e-17 1.20589352375036e-13
#MSTRG.15005 8.15398099367908 3.52136235857403e-16 1.92177691104112e-12

出现一MSTRG 开头的基因编号，它主要是因爲我之前在用hisat2做mapping 時用的基因組注釋文件和基因組index用的genome不一致導致的，我用自己下載的基因組做了一邊index後問題就解決了。

summary(res)
#out of 46687 with nonzero total read count
#adjusted p-value < 0.1
#LFC > 0 (up)       : 169, 0.36%
#LFC < 0 (down)     : 20, 0.043%
#outliers [1]       : 3020, 6.5%
#low counts [2]     : 7139, 15%
#(mean count < 5)
#[1] see 'cooksCutoff' argument of ?results
#[2] see 'independentFiltering' argument of ?results
##说明上述有169个基因上调，20个基因下调
 write.csv(res,file="All_results.csv")
 table(res$padj<0.05)
#FALSE  TRUE
#36366   162
##说明有padj<0.05的有162个。

提取差异表达genes(DEGs)并进行gene symbol注释

获取padj<0.05,表达倍数取log2的对数后>1或者<-1的差异表达基因。

diff_gene_deseq2 <- subset(res,padj < 0.05 & abs(log2FoldChange)>1)
dim(diff_gene_deseq2)
head(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file="DEG_treat_vs_control.csv")

6.2对差异表达基因进行注释


#(resOrdered <- res[order(res$padj), ])

#browseVignettes("DESeq2")

GO,KEGG,GSEA参数

7.GO/KEGG分析及GSEA(参考地址) 7.1安装Y 叔的R包clusterProfiler包，gitHub clusterProfiler是一个神器的包 Y叔的 clusterProfiler包的引用

G Yu, LG Wang, Y Han, QY He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology 2012, 16(5):284-287. doi:[10.1089/omi.2011.0118](http://dx.doi.org/10.1089/omi.2011.0118)

source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite("BiocUpgrade")
library(clusterProfiler)

biocLite("DOSE")
require(DOSE)
library(DO.db)

7.2安装构建自己的基因组注释数据 Biocouductor官网有19个物种，例如:载入人类的数据

biocLite("org.Hs.eg.db")
library(org.Hs.eg.db)

但是木有我要的。自己AnnotationHub构建自己的OrgDb。很巧合，Y叔写的有我的物种的一篇推送传送门另一篇参考地址。

必须提的准确性问题。徐州更提到的AnnotationHub的物种注释包的可靠性问题详细地址

library(AnnotationHub)
hub <- AnnotationHUb()

query(hub,"Zea mays")

s1 <- hub[[""]]

7.3 GO分析 [Y叔官方文档]((https://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/) ID转换，常见的ENSEMBL，ENTREZID。ENTREZID=kegg=ncbi-geneid ID转换函数 keytypes() gene <- row.names() 7.4 GSEA分析

基因集富集分析 (Gene Set Enrichment Analysis, GSEA) 的基本思想是使用预定义的基因集（通常来自功能注释或先前实验的结果），将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。基因集合富集分析检测基因集合而不是单个基因的表达变化，因此可以包含这些细微的表达变化，预期得到更为理想的结果。参考资料：GSEA分析是什么鬼(上)和GSEA分析是什么鬼(下)。

作者：lxmic 链接：https://www.jianshu.com/p/4910d7cec5c8

#Gene Set Enrichment Analysis (GSEA)

7.5 KEGG(pathway)分析

KEGG将基因组信息和高一级的功能信息有机地结合起来，通过对细胞内已知生物学过程的计算机化处理和将现有的基因功能解释标准化，对基因的功能进行系统化的分析。 KEGG的另一个任务是一个将基因组中的一系列基因用一个细胞内的分子相互作用的网络连接起来的过程，如一个通路或是一个复合物，通过它们来展现更高一级的生物学功能。参考文章：http://blog.sciencenet.cn/blog-364884-779116.html KEGG物种缩写：http://www.genome.jp/kegg/catalog/org_list.html GO和KEGG输出文件解读：http://www.bio-info-trainee.com/370.html 参考地址同上。

上述是在基因水平上的差异 下面找的是在转录本水平上的差异

Cufflinks

整个实验的大概脚本情况

#!/bin/bash
#批量循环，创造出文件名，并且复制文件。
#mkdir 102
A="CR"
B="Cr"
C="Nr"
tit[0]="_L8_I313."  #CR-1
tit[1]="L8_I314."   #CR-2
tit[2]="_L2_I301."  #Cr-1
tit[3]="L2_I302."           #Cr-2
tit[4]="_L2_I307."    #Nr-1
tit[5]="L2_I308."   #Nr-2

mid="-paired-"
end=".fastq.gz"
R1="-R1.fastq.gz"
R2="-R2.fastq.gz"

bam1="map.bam"
for (( i = 1; i < 4; i++ )); 
do
    k=$i
    # echo ${A}${mid}${i}${R1}
    # echo ${A}${mid}${i}${R2}
    # echo ${B}${mid}${i}${R1}
    # echo ${B}${mid}${i}${R2}
    # echo ${C}${mid}${i}${R1}
    # echo ${C}${mid}${i}${R2}

    #第1步质控：trimmomatic（共9条，更改输入输出文件名即可）
    #java -jar /home/guo/tool/Trimmomatic-0.38/trimmomatic-0.38.jar PE -threads 8 -phred33 Cr-zaoqian-1_L2_I301.R1.clean.fastq.gz Cr-zaoqian-1_L2_I301.R2.clean.fastq.gz Cr-paired-1-R1.fastq.gz Cr-unpaired-1-R1.fastq.gz Cr-paired-1-R2.fastq.gz Cr-unpaired-1-R2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 &

    #第2步比对：hisat2(共三条)
    #hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 ${A}${mid}${i}${R1} -2 ${A}${mid}${i}${R2} -S ${A}${i}".map.sam" --dta-cufflinks --novel-splicesite-outfile ${A}${i}".nsplice" &

    #hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 ${B}${mid}${i}${R1} -2 ${B}${mid}${i}${R2} -S ${B}${i}".map.sam" --dta-cufflinks --novel-splicesite-outfile ${B}${i}".nsplice" &

    #hisat2 -x /home/guo/maize/zm437/Zea_mays.AGPv4.dna.toplevel -p 8 -1 ${C}${mid}${i}${R1} -2 ${C}${mid}${i}${R2} -S ${C}${i}".map.sam" --dta-cufflinks --novel-splicesite-outfile ${C}${i}".nsplice" &
    #第3步：用samtool，格式转换，将sam转换为bam（共三条）
    # samtools sort -@ 8 -o ${A}${i}".map.bam" ${A}${i}".map.sam" &
    # samtools sort -@ 8 -o ${B}${i}".map.bam" ${B}${i}".map.sam" &
    # samtools sort -@ 8 -o ${C}${i}".map.bam" ${C}${i}".map.sam" &
    #第4步装配：用stringtie（共三轮）
    #第一轮（9个分别比对到基因组）
    #stringtie ${A}${i}".map.bam" -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o ${A}${i}".gtf" &
    #stringtie ${B}${i}".map.bam" -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o ${B}${i}".gtf" &
    #stringtie ${C}${i}".map.bam" -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o ${C}${i}".gtf" &
    #第二轮（整合9个的结果成一个）
    #stringtie --merge -G /home/guo/maize/zm437/Zea_mays.AGPv4.37.gtf -p 8 -o merged.gtf CR1.gtf CR2.gtf CR3.gtf Cr1.gtf Cr2.gtf Cr3.gtf Nr1.gtf Nr2.gtf Nr3.gtf &
    #第三轮（以第二轮的结果作为参考序列，9个分别比对）
    # mkdir ${A}${i}"_out"
    # mkdir ${B}${i}"_out"
    # mkdir ${C}${i}"_out"
    # stringtie ${A}${i}".map.bam" -G merged.gtf -p 8 -b ${A}${i}"_out" -e -o ${A}${i}"-st.gtf" &
    # stringtie ${B}${i}".map.bam" -G merged.gtf -p 8 -b ${B}${i}"_out" -e -o ${B}${i}"-st.gtf" &
    # stringtie ${C}${i}".map.bam" -G merged.gtf -p 8 -b ${C}${i}"_out" -e -o ${C}${i}"-st.gtf" &
    第5步 生成CSV文件
    #python路径 python2.7 /disks/backup/chaim/soft/prepDE.py -i 
    第6步 deseq2进行定量分析




















done