case2 - nibb-gitc/gitc2023mar-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
  • gene annotation: genes_except_chr5.gtf (genes.gtfからchr5上の情報を除いたもの)

Software

  • hisat2 (installed)
  • stringtie (installed)
  • samtools (installed)
  • edgeR (installed)
  • multiqc (installed)

Setup

Setup environment

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

$ mkdir case2
$ cd case2

Sequence reads

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

qualityを確認するなら以下も有効な情報が得られる。 seqkit stats -a *.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実行時にmappning stats情報がreportされる。これは記録しておくとよい。

あるいは以下、--new-summary --summary-fileのオプションを付けておくと、mappning stats情報が指定するファイルに出力される。 また後述のmultiqcで、さらにまとめたreportを得ることもできる。

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

Hisat2 stats summarize (参考）

multiqcを利用すると以下が可能である。

$multiqc -d .

で、カレントディレクトリーにあるstats情報を認識し、集計したhtmlファイル等が得られる。 https://multiqc.info/docs/#hisat2

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 で可視化。

  1. 手元のマシンに全ての sorted.bam を scp コマンドで転送してくること。

  2. IGVを立上げる。

  3. メニュー Genomes > Load Genome From File... で ｇenome.fa を選択、File > Load from File ... で genes.gtf を選択

  4. File > Load from File ... で作製した .sorted.bam を読み込む

  5. 適当にズームアップする。

Transcript assembly and quantification with StringTie

ここでは-eは付けず、-Gで指定したgtfファイルに記載の遺伝子モデルと共に、新規な遺伝子の発見も目指す。

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

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

出力されたgtfファイルのリストを作成 sample.list

stringtie_count_no_e/2D_rep1/2D_rep1.gtf
stringtie_count_no_e/2D_rep2/2D_rep2.gtf
stringtie_count_no_e/2D_rep3/2D_rep3.gtf
stringtie_count_no_e/2D2L_rep1/2D2L_rep1.gtf
stringtie_count_no_e/2D2L_rep2/2D2L_rep2.gtf
stringtie_count_no_e/2D2L_rep3/2D2L_rep3.gtf

ここではcommand user interfaceで以下の流れで作製

emacs sample.list
stringtie_count_no_e/2D_rep1/2D_rep1.gtf
stringtie_count_no_e/2D_rep2/2D_rep2.gtf
stringtie_count_no_e/2D_rep3/2D_rep3.gtf
stringtie_count_no_e/2D2L_rep1/2D2L_rep1.gtf
stringtie_count_no_e/2D2L_rep2/2D2L_rep2.gtf
stringtie_count_no_e/2D2L_rep3/2D2L_rep3.gtf
C-x C-c　←ctrlキーを押しながらxおよびc
y ←で確認

元のgene.gtfファイルに新規に見出されたgtfがmergeされたgtfファイルが出力される。

stringtie --merge -p 4 -G ../genes.gtf -o stringtie_merged_no_e.gtf sample.list

Inspect new anntation

gffcompareを使って、stringtie mergeによって追加されたgtf情報を確認

gffcompare -r ../genes.gtf -o merged stringtie_merged_no_e.gtf

merged.statsを見る、以下新規に229のexonが見出されている。

# gffcompare v0.12.6 | Command line was:
#gffcompare -r ../genes.gtf -o merged stringtie_merged_no_e.gtf
#

#= Summary for dataset: stringtie_merged_no_e.gtf 
#     Query mRNAs :   44126 in   33311 loci  (32630 multi-exon transcripts)
#            (7113 multi-transcript loci, ~1.3 transcripts per locus)
# Reference mRNAs :   41611 in   33350 loci  (30130 multi-exon)
# Super-loci w/ reference transcripts:    33274
#-----------------| Sensitivity | Precision  |
        Base level:   100.0     |    99.5    |
        Exon level:    99.6     |    97.5    |
      Intron level:   100.0     |    99.1    |
Intron chain level:   100.0     |    92.3    |
  Transcript level:    99.9     |    94.2    |
       Locus level:    99.9     |    99.8    |

     Matching intron chains:   30127
       Matching transcripts:   41551
              Matching loci:   33321

          Missed exons:       0/169268  (  0.0%)
           Novel exons:     229/171803  (  0.1%)
        Missed introns:       0/127896  (  0.0%)
         Novel introns:     232/129089  (  0.2%)
           Missed loci:       0/33350   (  0.0%)
            Novel loci:      37/33311   (  0.1%)

 Total union super-loci across all input datasets: 33311 
44126 out of 44126 consensus transcripts written in merged.annotated.gtf (0 discarded as redundant)

Transcript quantification with StringTie

再度、この新しいgtfファイルをrefとして、今度は-eを付加してカウントさせる。

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

Build count matrix

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

$ prepDE.py3 -i stringtie_count_e

Results => transcript_count_matrix.csv, gene_count_matrix.csv

less コマンドで確認

DEG test using edgeR

> library(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
MSTRG.1|ARV1                 8          8          4        5        0       19
AT1G01060|LHY                3          3          0        0        0        4
AT1G01070|AT1G01070         11          9          4        0        0        0
MSTRG.6|DCL1                36         18         22       24       42       18
MSTRG.6|MIR838A              0          0          1        0        0        0
MSTRG.7|AtPPa1              67         70         54       45       33       25

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
MSTRG.6096|LHB1B1   10.403846 10.606981 1.797278e-71 6.416821e-67
MSTRG.16054|PORA   -12.884845  9.322153 2.861336e-65 5.107913e-61
MSTRG.4846|AGT1     13.052784  9.491300 2.999409e-63 3.569597e-59
MSTRG.1971|AB165     8.280229 10.313082 7.898897e-62 7.050358e-58
MSTRG.16091|THI1     6.920644 10.675690 6.906838e-57 4.931897e-53
MSTRG.1364|ATLFNR2  12.134788  8.578639 6.279815e-56 3.294726e-52
MSTRG.8665|ICL      -6.486209 11.514574 6.459705e-56 3.294726e-52
MSTRG.11196|XTR7    -6.249156  9.324708 4.418952e-53 1.972123e-49
MSTRG.6393|LTP2     -5.451952 11.454657 5.350753e-53 2.122644e-49
MSTRG.9927|LHCA1     5.237290 11.382122 2.425320e-50 8.659121e-47

> 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)

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