RNA Seq 解析

MiSeqデータのリファレンス配列へのマッピングと転写産物量の推定。

3年ほど前に行っていた古い方法だが、一応アップしておきます。
これからやる方は、マッピングにはHISAT2を使ってください。

注：javaがインストールされていること。

クォリティーチェック（FastQC）
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

FastQCのインストールはホームページかららダウンロードして使うか、あるいはコマンドラインからダウンロードします。

$ wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.8.zip
$ unzip  fastqc_v0.11.2.zip   

FastQCフォルダの中のfastqcファイルを実行ファイルに変更

$ cd FastQC
$ chmod +x fastqc

圧縮されたシーケンスファイル（.fastaq.gz）ファイルの解凍
1) MacOSのFinder上でダブルクリック
2) gunzip
3) tar zxfv

クォリティーチェック（fastqc）は解凍してなくても可です。
例）20160801フォルダの中の1B2_S2_L001_R1_001.fastqファイルのクオリティーをチェックして、結果をFastQC_data/20160801_dataフォルダに出力する。

$ FastQC/fastqc --nogroup -o ./FastQC_data/20160801_data 20160801/1B2_S2_L001_R1_001.fastq

トリムは、fastq_quality_trimmer、Trimmomaticのどちらを使ってもよい。
（ペアエンドの場合はTrimmomaticの方がいいかも）

FASTX-Toolkitの中にあるfastq_quality_trimmerを利用する場合。
http://hannonlab.cshl.edu/fastx_toolkit/

$ fastq_quality_trimmer -Q 33 -t 20 -l 30 -i /Users/hanano/data/20160801/1B2_S2_L001_R1_001.fastq -o /Users/hanano/data/Trim_data/20160801_data/1B2_S2_L001_R1_001_trim.fastq

Trimmomaticを利用する場合。
http://www.usadellab.org/cms/?page=trimmomatic

シングルエンド

java -jar trimmomatic-0.33.jar \ SE \ # single-end -threads 4 \ # スレッド数 -phred33 \ # phred33 または -phred64 を指定 -trimlog log.txt \ # 実行ログの保存先 input.fq \ # 入力 FASTQ output.fq \ # 出力 FASTQ ILLUMINACLIP:adapters.fa:2:10:10 \ # アダプター除去条件の指定 MINLEN: 30 # 30bp を満たさないリードを除去

ペアエンド
java -jar trimmomatic-0.33.jar \ PE \ -threads 4 \ -phred33 \ -trimlog log.txt \ input_1.fq \ # 1つ目の FASTQ input_2.fq \ # 2つ目の FASTQ paired_output_1.fq \ # 1つ目の FASTQ のアダプター配列除去結果（paired output） unpaired_output_1.fq \ # 1つ目の FASTQ のアダプター配列除去結果（unpaired output） paired_output_2.fq \ # 2つ目の FASTQ のアダプター配列除去結果（paired output） unpaired_output_2.fq \ # 2つ目の FASTQ のアダプター配列除去結果（unpaired output） ILLUMINACLIP:adapters.fa:2:30:10 \ LEADING:20 \ TRAILING:20 \ SLIDINGWINDOW:4:15 \ MINLEN:36
java -jar trimmomatic-0.36.jar SE -phred33 sample_1B2R1.fastq sample_1B2R1_trim.fastq ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MIN
LEN:36

$ java -jar Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 2N6-4_S6_L001_R1_001.fastq.gz 2N6-4_S6_L001_R2_001.fastq.gz /Trimmomatic_results/2N6-4_S6_L001_R1_paired.fq.gz /Trimmomatic_results/2N6-4_S6_L001_R1_unpaired.fq.gz /Trimmomatic_results/2N6-4_S6_L001_R2_paired.fq.gz /Trimmomatic_results/2N6-4_S6_L001_R2_unpaired.fq.gz ILLUMINACLIP:Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

マッピング（tophat）　
以前はtophatで行なっていたが、現在tophatはメンテナンスが終了しているので、これから解析する場合にはHISAT2を使用してください。
http://ccb.jhu.edu/software/hisat2/index.shtml

$ tophat -p 2 -G Arabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Annotation/Archives/archive-2015-07-17-14-28-46/Genes/genes.gtf -o tophat_results/1B4_S6_L001_R1_001_trim_P0 Arabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Sequence/Bowtie2Index/genome 1B4_S6_L001_R1_001_trim.fastq

ペアエンドのファイルを一つの結果にマッピングする場合

tophat -p 4 -r 100 -G genes.gtf -o tophat_results bowtie_index_genome reads1_1.fastq reads1_2.fastq

$ tools/tophat-2.1.1.Linux_x86_64/tophat -p 2 -G Arabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Annotation/Archives/archive-2015-07-17-14-28-46/Genes/genes.gtf -o analysis/20170411/Tophat_results/2N6-4_S6_trim_P1 Arabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Sequence/Bowtie2Index/genome analysis/20170411/Trimmomatic_results/2N6-4_S6_L001_R1_paired.fq.gz analysis/20170411/Trimmomatic_results/2N6-4_S6_L001_R2_paired.fq.gz

遺伝子発現量の計算（cufflinks）
http://cole-trapnell-lab.github.io/cufflinks/

$ cufflinks -p 2 -g Arabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Annotation/Archives/archive-2015-07-17-14-28-46/Genes/genes.gtf tophat_results/1B2_S2_L001_R1_001_trim_P0/1B2_S2_L001_R1_001_trim_P0.bam -o tophat_results/1B2_S2_L001_R1_001_trim_P0/cufflinks_results

$ ls analysis/20180516/tophat_results/2A-0_S2_trim/cufflinks_results/2A-0_S2_transcripts.gtf analysis/20180516/tophat_results/2A-D3_S6_trim/cufflinks_results/2A-D3_S6_transcripts.gtf analysis/20180516/tophat_results/2A-E3_S7_trim/cufflinks_results/2A-E3_S7_transcripts.gtf analysis/20180516/tophat_results/C-D1_S3_trim/cufflinks_results/C-D1_S3_transcripts.gtf analysis/20180516/tophat_results/C-E1_S4_trim/cufflinks_results/C-E1_S4_transcripts.gtf analysis/20180516/tophat_results/2A-D1_S5_trim/cufflinks_results/2A-D1_S5_transcripts.gtf analysis/20180516/tophat_results/2A-E1_S6_trim/cufflinks_results/2A-E1_S6_transcripts.gtf analysis/20180516/tophat_results/C-0_S1_trim/cufflinks_results/C-0_S1_transcripts.gtf analysis/20180516/tophat_results/C-D2_S7_trim/cufflinks_results/C-D2_S7_transcripts.gtf analysis/20180516/tophat_results/C-E2_S1_trim/cufflinks_results/C-E2_S1_transcripts.gtf analysis/20180516/tophat_results/2A-D2_S2_trim/cufflinks_results/2A-D2_S2_transcripts.gtf analysis/20180516/tophat_results/2A-E2_S3_trim/cufflinks_results/2A-E2_S3_transcripts.gtf analysis/20180516/tophat_results/C-D3_S4_trim/cufflinks_results/C-D3_S4_transcripts.gtf analysis/20180516/tophat_results/C-E3_S5_trim/cufflinks_results/C-E3_S5_transcripts.gtf > transcripts.gtf.txt

結果の統合（cuffmerge)

$ tools/cufflinks-2.2.1.Linux_x86_64/cuffmerge  -p 2 -o analysis/20170411/fpkm_compare/2C3-2 -g Arabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Annotation/Archives/arArabidopsis_thaliana_NCBI_TAIR10/Arabidopsis_thaliana/NCBI/TAIR10/Sequence/Bowtie2Index/genome.fa analysis/20170411/2C3-2_transcripts.gtf.txt

cuffdiff

$ /home/hanano/tools/cufflinks-2.2.1.Linux_x86_64/cuffdiff -p 2 --upper-quartile-norm -o /home/hanano/analysis/20170411/cuffdiff_result -L 2C3-2,2C3-3,2C3-4,2C6-2,2C6-3,2C6-4,2N3-2,2N3-3,2N3-4,2N6-2,2N6-3,2N6-4 /home/hanano/analysis/20170411/fpkm_compare/merged.gtf /home/hanano/analysis/20170411/Tophat_results/2C3-2_S1_trim_P1/2C3-2_S1_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2C3-3_S5_trim_P1/2C3-3_S5_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2C3-4_S3_trim_P1/2C3-4_S3_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2C6-2_S3_trim_P1/2C6-2_S3_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2C6-3_S1_trim_P1/2C6-3_S1_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2C6-4_S5_trim_P1/2C6-4_S5_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2N3-2_S2_trim_P1/2N3-2_S2_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2N3-3_S6_trim_P1/2N3-3_S6_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2N3-4_S4_trim_P1/2N3-4_S4_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2N6-2_S4_trim_P1/2N6-2_S4_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2N6-3_S2_trim_P1/2N6-3_S2_trim_P1.bam /home/hanano/analysis/20170411/Tophat_results/2N6-4_S6_trim_P1/2N6-4_S6_trim_P1.bam

* トリム以降は、時間がかかるのでshellスクリプトを書いてqsubでスパコンに投げた方がよい。
qsubを使ったjobの投入は、基本的には’**.sh’ファイルを作成して、’qsub sh **.sh’だが、メモリの設定などオプションが様々なので、管理者の指定に従うこと。

例）

$ qsub -l s_vmem=6G,mem_req=6G -soft -l short,month -cwd -S /bin/sh Trimmomatic_20180205.sh
 
$ qsub -jc k_P1 -S /bin/bash -e error.log -cwd -b y -o qsub_out sh cuffdiff_20170502.sh