番茄轉錄組高通量數據分析學習筆記 - Part I

01-30

前言

Hello，大家好！

我是米粒兒，我是一名做園藝類實驗的研究生，這是我第一次做高通量的數據分析。我大概利用1個月不到的時間完成了第一次RNA-Seq數據分析流程的跑通，這是我的學習筆記的第一部分。至於第二部分的GO分析，KEGG Pathway分析等等我還在繼續整理，希望大家能夠喜歡！

另外，我在本文的結尾給大家推薦一下講生物信息學分析的知乎Live課程，我的入門是完全通過聽課完成的（這是硬廣）。

練習數據來源文獻

Transcriptome Analysis of the Cf-12-Mediated Resistance Response to Cladosporium fulvum in Tomato，Front. Plant Sci., 05 January 2017 | (https://doi.org/10.3389/fpls.2016.02012)

數據背景：三個時間0dpi,4dpi和8dpi代表三個處理，每個處理三個rep

下載數據

文獻中提供的數據信息：The raw sequencing data of the nine samples have been submitted to the NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra). The accession numbers are: SRR4041970, SRR4041973, SRR4041974, SRR4041975, SRR4042017, SRR4042029, SRR4042030, SRR4042031, and SRR4042032.

下載數據使用的命令:

~$ prefetch –v SRR4041970n

ps. 下載 SRR4041970.sra 的文件默認安裝到 sratoolkit 時配置的 public 目錄中，位置是 ~/ncbi/public/sra

重命名sra文件:

mv SRR4041970.sra Cf-12-RNA-Seq_0dpi_rep1_paired.sran

ps.同文件夾內改名不需要備註路徑。

將sra文件改為fastqc.gz格式修改文件格式

fastq-dump --split-files --gzip ./raw.sra/Cf-12-RNA-Seq_0dpi_rep1_paired.sra --outdir ./raw.fastq &n

ps. fastq-dump, 將sra文件更改為fastq.gz 格式；--gzip, 輸出文件壓縮成gzip格式; --split-files，針對pair-end測序結果，將原文件轉換為兩個fastq.gz文件，裡面的reads是成對的；&, 後台運行;

對fastq.gz文件進行fastqc,檢測reads質量

fastqc -q -t 18 ./raw.fastq/*.fastq.gz -o ./fastqc.raw &n

ps. fastqc, 評價測序質量，輸出結果到fastqc.raw；-q,安靜運行模式；-t，選擇程序運行線程數；-o，fastqc生成的報告文件的儲存路徑.

fastqc結果分析

具體參考孟大神 https://zhuanlan.zhihu.com/p/20731723 非常詳細

根據fastqc結果切除cutadaptor

根據上一步fastqc結果中adaptor content結果，使用cutadaptor對reads進行處理

nohup cutadapt --times 1 -e 0.1 -O 3 -m 115 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o $/Cf-12-RNA-Seq_0dpi_rep1_paired_1_cut.fastq.gz -p $/Cf-12-RNA-Seq_0dpi_rep1_paired_2_cut.fastq.gz $/Cf-12-RNA-Seq_0dpi_rep1_paired_1.fastq.gz $/Cf-12-RNA-Seq_0dpi_rep1_paired_2.fastq.gz > $/Cf-12-RNA-Seq_0dpi_rep1_peired_cutadapt.log 2>&1 &n

PS, --times,切除次數；-e 0.1,match時允許10%的錯誤率；-o 3，序列與adaptor序列不少於3個base的match；-m,cut adaptor之後當序列長度小於115的時候該序列捨棄；-o,輸出的reads1的文件；-p，輸出的reads2的文件

切除完adaptor之後，根據fastqc結果切除reads首尾測序質量較差的片段

一般我們確定要trim多少bp的時候都是按phred20這個標準進行評估，調用fastx_trimmer

zcat ./raw.fastq/Cf-12-RNA-Seq_0dpi_rep1_paired_1_cut.fastq.gz | fastx_trimmer -f 16 -l 115 -z -o ./raw.fastq/Cf-12-RNA-Seq_0dpi_rep1_paired_1_cut_trimmer.fastq.gz &n

ps, zcat，將gz文件解壓縮；"|",解壓的文件被fastx_trimmer調用；-f,保留的第一個base，-l,保留的最後一個base；-z,將輸出文件進行壓縮；

使用bowtie2分析reads中mapping到rRNA的比例，並使用tophat2將unmap數據比對到參考基因組上

在RNA-Seq分析流程中，使用bowtie2或者tophat2將reads比對到參考rRNA上能夠分析建庫過程中是否存在rRNA的混入，真核生物rRNA非常保守，所以其參考序列是一定的（番茄的rRNA包括25S，17S和4.5S)，需要用bowtie-build building a new index

NCBI下載; 搜索 "rRNA [Feature key]tomato",下載tomato fasta格式17S,25S,4.5S的rRNA序列,由transmit上傳到伺服器；
merge fasta 文件 to tomato_rRNA.fa

cat ./*.fasta > ./tomato_rRNA.fan

build tomato_rRNA_index

boetie2-build ./tomato_rRNA.fa ./tomato_rRNAn

結果

tomato_rRNA.1.bt2
tomato_rRNA.2.bt2
tomato_rRNA.3.bt2
tomato_rRNA.4.bt2
tomato_rRNA.rev.1.bt2
tomato_rRNA.rev.2.bt2
調用tomato_rRNA

利用bowtie2對reads進行比對

nohup bowtie2 -p 5 -S $/Cf-12-RNA-Seq_0dpi_rep1_bowtie -x $/tomato_rRNA -1 $/Cf-12-RNA-Seq_0dpi_rep1_paired_1_cut_trimmer.fastq.gz -2 $/Cf-12-RNA-Seq_0dpi_rep1_paired_2_cut_trimmer.fastq.gz --un-conc $/Cf-12-RNA-Seq_0dpi_rep1_unmap > $/Cf-12-RNA-Seq_0dpi_rep1_bowtie.log 2>&1 &n

ps, -p,調用核心數目；-S，輸出文件到；-x,bt2-idx; -1,paired_1; -2, paired_2; --un-conc,輸出unmap的文件到；
三種格式，.fasta = .fa = .fsa

下載番茄參考基因組文件並使用bowtie2-bulid新建genome_index

下載植物參考基因組推薦網站

phytozome https://phytozome.jgi.doe.gov/pz/portal.html
silva https://www.arb-silva.de/

使用tophat2將unmap數據比對到參考基因組上

nohup tophat2 -p 5 -o ./tophat2-result/Cf-12-RNA-Seq_0dpi_rep1 $/SL2.50_genome $/bowtie2_rRNA_result/Cf-12-RNA-Seq_0dpi_rep1_unmap.1 $/bowtie2_rRNA_result/Cf-12-RNA-Seq_0dpi_rep1_unmap.2 > ./tophat2-result/Cf-12-RNA-Seq_0dpi_rep1.log 2>&1 &n

使用cufflinks拼接轉錄本並估算表達丰度

nohup cufflinks -o ./cufflink_result/Cf-12-RNA-Seq_0dpi_rep1_cufflink -p 3 -G /home/zhangyt/reference/SL2.50_gff/Solanum_lycopersicum.SL2.50.37.chr.change_name.gff3 ./tophat2-result/Cf-12-RNA-Seq_0dpi_rep1/accepted_hits.bam > ./cufflink_result/Cf-12-RNA-Seq_0dpi_rep1_cufflink.log 2>&1 &n

ps,-o,輸出文件到；-G，需要用的轉錄組參考文件；

使用cuffdiff計算比較不同處理間基因表達差異(以0dpi和4dpi的處理比較為實例)

nohup cuffdiff -o ./cuffdiff_result/Cf-12-RNA-Seq_0dpi_4dpi_cuffdiff -p 16 --labels Cf-12_0dpi,Cf-12_4dpi --min-reps-for-js-test 3 $/SL2.50_gff/Solanum_lycopersicum.SL2.50.37.chr.change_name.gff3 ${Cf_0dpi_bam} ${Cf_4dpi_bam} > ./cuffdiff_result/Cf-12-RNA-Seq_0dpi_4dpi_cuffdiff.log 2>&1 &n

ps, -p,d調用核心數目；--lables,輸入的文件次序（先0後4）;--min-reps-for-js-test，每個處理的rep數目；.gff3,番茄基因組的注釋文件；

Linux語言小注釋

wget 下載文件
mkdir 新建文件夾
vim，編輯shell文件； bash，運行shell文件
linux 中的文件名不能出現空格，否則在後續使用過程中會出現問題