case2 - nibb-gitc/gitc2021mar-rnaseq GitHub Wiki

Case study 2: Genome-based RNA-Seq pipeline

アラビドプシス(Arabidopsis thaliana)のRNA-seqを行った。ライブラリは2D sample (2days dark conditionで生育させた黄色芽生え)と2D2L sample (その後さらに2days light conditionで生育させた緑化芽生え）でそれぞれsampling duplicateを３つ用意した。シーケンスはpaired-end（インサートの両端を読む）で76bpシークエンスしたものを事前にpre-processing (cutadapt処理済)している。

これらのdataを用いて、HISAT2 -> StringTie -> edgeR のパイプラインで DEG (differential expressed genes)を明らかにせよ。

マッピング結果をIGVで可視化し、DEGのを確認せよ。

Data

Input reads

(ファイルは、~/gitc/data/KY/genome_base/　にある)

condition Dark, rep#1: 2D_rep1_R1.fastq, 2D_rep1_R2.fastq [R1: Read1 (Fwd); R2: Read2 (Rev) 以下同様]
condition Dark, rep#2: 2D_rep2_R1.fastq, 2D_rep2_R2.fastq
condition Dark, rep#3: 2D_rep3_R1.fastq, 2D_rep3_R2.fastq
condition Light, rep#1: 2D2L_rep1_R1.fastq, 2D2L_rep1_R2.fastq
condition Light, rep#2: 2D2L_rep2_R1.fastq, 2D2L_rep2_R2.fastq
condition Light, rep#3: 2D2L_rep3_R1.fastq, 2D2L_rep3_R2.fastq

Reference

genome sequence: genome.fa
gene annotation: genes.gtf

Software

hisat2 (installed)
stringtie (installed)
samtools (installed)
edgeR (installed)

Setup

Setup environment

case2 ディレクトリをつくり、以下の解析はその下で作業しよう。

$ mkdir case2
$ cd case2

Sequence reads

''less'' などのコマンドで、 2D_rep1_R1.fastq の内容を確認する。

注）Pre-processing済みであることが分かる。

Mapping with hisat2

Build hisat2 index from reference genome

hisat2-build genome.fa genome

Run hisat2

$ hisat2 -p 4 --dta -x genome -1 2D_rep1_R1.fastq -2 2D_rep1_R2.fastq -S 2D_rep1.sam

他の５つも同様に。

$ hisat2 -p 4 --dta -x genome -1 2D_rep2_R1.fastq -2 2D_rep2_R2.fastq -S 2D_rep2.sam
$ hisat2 -p 4 --dta -x genome -1 2D_rep3_R1.fastq -2 2D_rep3_R2.fastq -S 2D_rep3.sam
$ hisat2 -p 4 --dta -x genome -1 2D2L_rep1_R1.fastq -2 2D2L_rep1_R2.fastq -S 2D2L_rep1.sam
$ hisat2 -p 4 --dta -x genome -1 2D2L_rep2_R1.fastq -2 2D2L_rep2_R2.fastq -S 2D2L_rep2.sam
$ hisat2 -p 4 --dta -x genome -1 2D2L_rep3_R1.fastq -2 2D2L_rep3_R2.fastq -S 2D2L_rep3.sam

*-p は使うCPU coreを指定するオプション。使用するコンピュータのスペックに合わせて。

Inspect Results

hisat2実行時にmapping rateがレポートされる。これは記録しておくとよい。

Convert SAM to BAM

samtools sort -@ 4 -o 2D_rep1.sorted.bam  2D_rep1.sam
samtools index 2D_rep1.sorted.bam

他の５つも同様。ちなみに、シェルスクリプトで６つ自動で処理するには以下のようにする。

for f in *sam
do 
  samtools sort -o `basename $f .sam`.sorted.bam $f
done

for f in *bam
do
  samtools index $f
done

IGV

IGV で可視化。

手元のマシンに全ての sorted.bam を scp コマンドで転送してくること。
IGVを立上げる。
メニュー Genomes > Load Genome From File... で ｇenome.fa を選択、File > Load from File ... で genes.gtf を選択
File > Load from File ... で作製した .sorted.bam を読み込む
適当にズームアップする。

Counting with StringTie

$ stringtie -e -p 4 -G genes.gtf -o stringtie_count/2D_rep1/2D_rep1.gtf 2D_rep1.sorted.bam
$ stringtie -e -p 4 -G genes.gtf -o stringtie_count/2D_rep2/2D_rep2.gtf 2D_rep2.sorted.bam
$ stringtie -e -p 4 -G genes.gtf -o stringtie_count/2D_rep3/2D_rep3.gtf 2D_rep3.sorted.bam
$ stringtie -e -p 4 -G genes.gtf -o stringtie_count/2D2L_rep1/2D2L_rep1.gtf 2D2L_rep1.sorted.bam
$ stringtie -e -p 4 -G genes.gtf -o stringtie_count/2D2L_rep2/2D2L_rep2.gtf 2D2L_rep2.sorted.bam
$ stringtie -e -p 4 -G genes.gtf -o stringtie_count/2D2L_rep3/2D2L_rep3.gtf 2D2L_rep3.sorted.bam

less コマンドstringtie_count/2D2L_rep2/2D2L_rep2.gtf 等gtfの中身を確認しておく。

Build count matrix

prepDE.py を使う。(stringtie開発者がedgeR / DEseq2用に提供)

$ python2 prepDE.py -i stringtie_count

Results => transcript_count_matrix.csv, gene_count_matrix.csv

less コマンドで確認

DEG test using edgeR

> dat <- read.csv("gene_count_matrix.csv", row.names=1)
> head(dat)
          X2D2L_rep1 X2D2L_rep2 X2D2L_rep3 X2D_rep1 X2D_rep2 X2D_rep3
AT4G22890        296        204        204       20       22       17
AT1G38440          0          0          0        0        0        0
AT3G27910          0          0          0        0        0        0
AT1G06620          3          0          6        0        3        4
AT5G54067          0          0          0        0        0        0
AT2G34630         52         13         10        9        0        3

group <- c(rep("2D2L", 3), rep("2D", 3))
D <- DGEList(dat, group=group)
D <- calcNormFactors(D)
D <- estimateCommonDisp(D)
D <- estimateTagwiseDisp(D)

> de.2D.vs.2D2L <- exactTest(D, pair=c("2D", "2D2L"))
> topTags(de.2D.vs.2D2L)
Comparison of groups:  2D2L-2D
               logFC    logCPM        PValue           FDR
AT2G34430  10.407862 10.708836 2.722476e-110 9.148063e-106
AT5G54190 -12.947382  9.462160 7.242338e-103  1.216785e-98
AT3G01500   4.707974 11.032867  9.737231e-98  1.090635e-93
AT2G38530  -5.485669 11.601129  1.773318e-97  1.489675e-93
AT5G54770   6.869227 10.780624  3.379061e-91  2.270864e-87
AT3G21720  -6.540733 11.652207  7.167374e-88  4.013968e-84
AT3G54890   5.097595 11.400061  2.520161e-86  1.209749e-82
AT1G73120 -11.504585  8.026029  3.565484e-80  1.497592e-76
AT4G14130  -6.562635  9.461145  3.723257e-79  1.390099e-75
AT2G05070   5.396477  9.093835  2.557098e-77  8.592362e-74

> tmp <- topTags(de.2D.vs.2D2L, n=nrow(de.2D.vs.2D2L))
> write.table(tmp$table, "de.2D.vs.2D2L.txt", sep="\t", quote=F)

発展問題

Q: How many genes are deferentially expressed?

Q: Scatter plot やMA plotを書いてみよう

Links

hisat2

https://ccb.jhu.edu/software/hisat2/index.shtml

stringtie

http://www.ccb.jhu.edu/software/stringtie/

edgeR

http://bioconductor.org/packages/release/bioc/html/edgeR.html

Notes

論文情報

Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown https://www.nature.com/articles/nprot.2016.095