Y叔的RNA-seq workflow 強力推薦
生信技能樹講解
- 數(shù)據(jù)已經(jīng)存在服務(wù)器，基因組注釋文件舶治、基因組文件都已經(jīng)存在箍邮。服務(wù)器環(huán)境軟件都已經(jīng)安裝完成昨悼。但是沒有root權(quán)限并扇。
數(shù)據(jù)獲取預(yù)處理
find ./510 -name '*zaoqian*'|xargs -i mv {} ./100 #批量查找510目錄里所有包含zaoqian的文件，移動到100目錄里。
- 實驗流程和使用軟件
質(zhì)控 trimmomatic
trimmomatic功能講解
精華篇 trimmomatic參數(shù)講解定量分析 Deseq2
轉(zhuǎn)錄本分析 Cufflink
從頭組裝 trinity
- 開始第一步——質(zhì)控。
fastqc可以檢測測序質(zhì)量血巍。
fastqc *.fastq.gz -t 4 #開啟4個進程。
trimmomatic去除接頭和低質(zhì)量測序數(shù)據(jù)（-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 &
最后的&是在后臺執(zhí)行珊随。
- 比對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之后禀崖，服務(wù)器實在吃不消衩辟。所以就不添加snp位點了。內(nèi)存在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 &
參數(shù)講解
-x /disks/backup/chaim/maize/genome_tran #基因組索引文件(后綴是.fai,此處只寫文件名艺晴，不要.fai后綴)
-p 8 #使用8個進程
-1 paired1 #質(zhì)控后左鏈的數(shù)據(jù)cms
-2 paired2 #質(zhì)控后右鏈的數(shù)據(jù)
-S Cr-1.map.sam #輸出文件
-dta-cufflinks--novel-splicesite-outfile Cr-1.nsplice #后期使用cufflinks做準備。
4.5 samtools將sam轉(zhuǎn)bam同時排序掸屡。
samtools sort -@ 8 -o ${bef[$i]}".map.bam" ${bef[$i]}".map.sam" 2>${bef[$i]}"samtool_out"
#構(gòu)建索引文件封寞，方便使用IGV查看拼接結(jié)果
samtools index -@ 8 {bef[$i]}".map.bam" {bef[$i]}".map.bai" &
- stringtie官網(wǎng)文檔（自帶梯子）
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
參數(shù)講解：
-G 參考基因組注釋文件 （我使用 gtf格式的注釋，到第三次會出錯仅财，所以使用的是gff3的注釋）
-p 8個進程
-o 輸出注釋文件名
-b 輸出其他文件的路徑（會同時生成多個文件狈究，一定要輸出在不同的路徑中，不然后面的輸出會覆蓋前面的輸出盏求。）
5.1 stringtie的其他功能——gffcompare和gffread
生信技能樹參考
官方地址
6.格式轉(zhuǎn)換抖锥，數(shù)據(jù)整合。
使用的是python的腳本碎罚。下載地址
python2.7 /disks/backup/chaim/soft/prepDE.py -i gtf2
#注意此處的prepDE.py的python版本是python2.7磅废，一定要指定版本號。否則會報錯荆烈。
gtf2文件內(nèi)容是對應(yīng)的比對后文件名和文件的位置拯勉。
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
輸出文件為
轉(zhuǎn)錄本輸出矩陣 transcript_count_matrix.csv
基因輸出矩陣 gene_count_matrix.csv
[perpDE.py參數(shù)](http://ccb.jhu.edu/software/stringtie/index.shtml?
t=manual#deseq)
-i 輸入文件
-g gene輸出矩陣位置
-t 轉(zhuǎn)錄本輸出矩陣位置
- Deseq2定量分析參考地址
有重復(fù)的差異分析使用DESeq2,
無重復(fù)的差異分析DEGseq.
BiocManager::install("DEGseq")
安裝到服務(wù)器端
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為基準比較上調(diào)或下調(diào)谜喊。deseq2默認的levels是c("未處理","處理")一定要注意這個順序潭兽。不然可能會產(chǎn)生相反的上升和下降值。
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
##上面構(gòu)建colData和countData,應(yīng)根據(jù)自己實際樣本數(shù)略作調(diào)整斗遏。但應(yīng)始終保持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
出現(xiàn)一MSTRG 開頭的基因編號诵次，它主要是因爲我之前在用hisat2做mapping 時用的基因組注釋文件和基因組index用的genome不一致導(dǎo)致的账蓉，我用自己下載的基因組做了一邊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個基因上調(diào)逾一，20個基因下調(diào)
write.csv(res,file="All_results.csv")
table(res$padj<0.05)
#FALSE TRUE
#36366 162
##說明有padj<0.05的有162個铸本。
提取差異表達genes(DEGs)并進行g(shù)ene symbol注釋
獲取padj<0.05,表達倍數(shù)取log2的對數(shù)后>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參數(shù)
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
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite("BiocUpgrade")
library(clusterProfiler)
biocLite("DOSE")
require(DOSE)
library(DO.db)
7.2安裝構(gòu)建自己的基因組注釋數(shù)據(jù)
Biocouductor官網(wǎng)有19個物種箱玷，
例如:載入人類的數(shù)據(jù)
biocLite("org.Hs.eg.db")
library(org.Hs.eg.db)
但是木有我要的。自己AnnotationHub構(gòu)建自己的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轉(zhuǎn)換舶得，常見的ENSEMBL，ENTREZID爽蝴。ENTREZID=kegg=ncbi-geneid
ID轉(zhuǎn)換函數(shù)
keytypes()
gene <- row.names()
7.4 GSEA分析
基因集富集分析 (Gene Set Enrichment Analysis, GSEA) 的基本思想是使用預(yù)定義的基因集（通常來自功能注釋或先前實驗的結(jié)果）沐批，將基因按照在兩類樣本中的差異表達程度排序，然后檢驗預(yù)先設(shè)定的基因集合是否在這個排序表的頂端或者底端富集蝎亚【藕ⅲ基因集合富集分析檢測基因集合而不是單個基因的表達變化，因此可以包含這些細微的表達變化颖对，預(yù)期得到更為理想的結(jié)果捻撑。
參考資料：GSEA分析是什么鬼(上)和GSEA分析是什么鬼(下)。
作者：lxmic
鏈接：http://www.reibang.com/p/4910d7cec5c8
#Gene Set Enrichment Analysis (GSEA)
7.5 KEGG(pathway)分析
KEGG將基因組信息和高一級的功能信息有機地結(jié)合起來缤底，通過對細胞內(nèi)已知生物學(xué)過程的計算機化處理和將現(xiàn)有的基因功能解釋標準化顾患，對基因的功能進行系統(tǒng)化的分析。
KEGG的另一個任務(wù)是一個將基因組中的一系列基因用一個細胞內(nèi)的分子相互作用的網(wǎng)絡(luò)連接起來的過程个唧，如一個通路或是一個復(fù)合物江解，通過它們來展現(xiàn)更高一級的生物學(xué)功能。
參考文章：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
參考地址同上徙歼。
上述是在基因水平上的差異
下面找的是在轉(zhuǎn)錄本水平上的差異
- Cufflinks
整個實驗的大概腳本情況
#!/bin/bash
#批量循環(huán)犁河，創(chuàng)造出文件名鳖枕，并且復(fù)制文件。
#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步質(zhì)控：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宾符，格式轉(zhuǎn)換，將sam轉(zhuǎn)換為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個的結(jié)果成一個）
#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 &
#第三輪（以第二輪的結(jié)果作為參考序列灭翔，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
