kallisto (T. aestivum, paired-end RNA-Seq) | kallisto を利用した 6 倍体小麦の転写産物発現量の定量

kallisto（Bray et al, 2016）を使って、6 倍体のパン小麦（T. aestivum）の paired-end RNA-Seq データ（Ma et al, 2017）を、Chinese Spring (IWGSC RefSeq v1.0) のリファレンスゲノム（IWGSC, 2018）にマッピングしてみる。ここで、乾燥ストレス実験で取られた RNA-Seq データを利用する。データはは、PRJNA380841 でアーカイブされている。

RNA-Seq データ

データベースにアーカイブされた PRJNA380841 のデータを利用する。データは全部で 24 サンプルからなるが、HISAT2 の実行テストであれば、数サンプルをダウンロードすればよい。以下の例では、bash を使って 24 サンプル 48 ファイルをローカルにダウンロードしている。

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR538/000/${seqlib}/${seqlib}_1.fastq.gz
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR538/000/${seqlib}/${seqlib}_2.fastq.gz
done

ダウンロードしたシーケンサーデータ（FASTQ ファイル）に対して、低クオリティリードの除去などのフィルタリングを行う。フィルタリング後のファイルを *_1.clean.fastq.gz および *_2.clean.fastq.gz のような名前で保存する。ここでは Trimmomatic を使用する。ファイルが多いので、bash の for 文で処理する。

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

for seqlib in ${SEQLIBS[@]}; do
    java -jar trimmomatic-0.38.jar PE -phred33 -threads 4    \
    ${seqlib}_1.fastq.gz ${seqlib}_2.fastq.gz                \
    ${seqlib}_1.clean.fastq.gz ${seqlib}_1.unpaired.fastq.gz \
    ${seqlib}_2.clean.fastq.gz ${seqlib}_2.unpaired.fastq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:60
done

ゲノムの配列データ

小麦（T. aestivum）Chinese Spring 系統のゲノム配列は URGI の FTP サイトで入手できる。wget コマンドで、ゲノム配列をダウンロードしてくる。ダウンロードした ZIP を展開すると、iwgsc_refseqv1.0_all_chromosomes フォルダが作られる。このフォルダの中に、小麦 Chinese Spring のゲノム配列が FASTA 形式で保存されている。

wget https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/iwgsc_refseqv1.0_all_chromosomes.zip
unzip iwgsc_refseqv1.0_all_chromosomes.zip
ls iwgsc_refseqv1.0_all_chromosomes
## 161010_Chinese_Spring_v1.0_pseudomolecules_AGP.tsv  161010_Chinese_Spring_v1.0_pseudomolecules_parts.fasta       161010_Chinese_Spring_v1.0_pseudomolecules_parts_to_scaffolds.bed  Agreement.txt
## 161010_Chinese_Spring_v1.0_pseudomolecules.fasta    161010_Chinese_Spring_v1.0_pseudomolecules_parts_to_chr.bed  161010_Chinese_Spring_v1.0_pseudomolecules_to_scaffolds.bed        README

kallisto の定量では、ゲノム配列を利用するのではなく、転写産物の配列を利用する。そのため、ゲノムのアセンブリ配列から転写産物の配列を抽出する必要がある。ここで、小麦のアノテーションファイルをダウンロードして、転写産物配列の抽出に利用する。現在、Chinese Spring のゲノム配列 RefSeq v1.0 に対して、アノテーションファイルは v1.0 と v1.1 の 2 種類提供されているが、ここでは新しいバージョンをダウンロードしておく。ダウンロードした ZIP を展開すると、アノテーションファイル（GFF/GTF）が展開されたフォルダに保存されている。

wget https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.1/iwgsc_refseqv1.1_genes_2017July06.zip
unzip iwgsc_refseqv1.1_genes_2017July06.zip
ls iwgsc_refseqv1.1_genes_2017July06
## IWGSC_v1.1_HC_20170706_cds.fasta  IWGSC_v1.1_HC_20170706.idmapping.csv      IWGSC_v1.1_LC_20170706_cds.fasta  IWGSC_v1.1_LC_20170706_pep.fasta          IWGSC_v1.1_LC_removed.ids
## IWGSC_v1.1_HC_20170706.gff3       IWGSC_v1.1_HC_20170706_pep.fasta          IWGSC_v1.1_LC_20170706.gff3       IWGSC_v1.1_LC_20170706_transcripts.fasta  README.txt
## IWGSC_v1.1_HC_20170706.gtf        IWGSC_v1.1_HC_20170706_transcripts.fasta  IWGSC_v1.1_LC_20170706.gtf        IWGSC_v1.1_LC.correspondanceTEHC.txt

次に Cufflinks 中の gffread コマンドを利用して、遺伝子アノテーション（IWGSC_v1.1_HC_20170706.gff3）を利用して、ゲノム配列（161010_Chinese_Spring_v1.0_pseudomolecules.fasta）から転写産物配列を取り出して、transcript.hc.fasta ファイルに保存する。

gffread -w transcript.hc.fasta -g 161010_Chinese_Spring_v1.0_pseudomolecules.fasta IWGSC_v1.1_HC_20170706.gff3

kallisto によるマッピング

kallisto でマッピングを行う前に、上で作成した転写産物の配列（transcript.hc.fasta）に対してインデックスを作成しておく。インデックスの名前を任意につけられるが、ここでは IGWSCv1_kallisto とした。

kallisto index -i IGWSCv1_kallisto transcript.hc.fasta

ls -h IGWSCv1_kallisto
## -rw-r--r-- 1 ____ ________ 2.7G Aug 27 11:32 IGWSCv1_kallisto

これが正しく実行されると、拡張子なしの IGWSCv1_kallisto ファイルが生成される。このファイルが kallisto が使用するインデックスとなる。

次に、kallisto を使用して転写産物の定量を行なっていく。定量結果の保存先フォルダは -o で指定する。ここでは SRR5383650_exp_kallisto のように保存先フォルダを指定した。-i に上で作成したインデックスを指定する。また、定量後に sleuth による発現変動遺伝子の検出を行う予定がある場合、-b オプションをつけて、ブートストラップ法による発現量の不確実性（「分散・ばらつき」のようなもの）を推定する。

SEQLIBS=(SRR5383650 SRR5383900 SRR5383901 SRR5383933 SRR5383934 SRR5383935 SRR5383936 SRR5383938 SRR5383941 SRR5383944 SRR5383945 SRR5383947 SRR5383948 SRR5383950 SRR5383952 SRR5383954 SRR5383955 SRR5383957 SRR5383959 SRR5383960 SRR5383961 SRR5383962 SRR5383963 SRR5383964)

for seqlib in ${SEQLIBS[@]}; do
    result_dir=${seqlib}_exp_kallisto
    kallisto quant -i IGWSCv1_kallisto -o ${result_dir} -b 100 ${seqlib}_1.clean.fastq.gz ${seqlib}_2.clean.fastq.gz 
done
##
## [quant] fragment length distribution will be estimated from the data
## [index] k-mer length: 31
## [index] number of targets: 133,744
## [index] number of k-mers: 121,817,028
## [index] number of equivalence classes: 588,686
## [quant] running in paired-end mode
## [quant] will process pair 1: SRR5383650_1.clean.fastq.gz
##                              SRR5383650_2.clean.fastq.gz
## [quant] finding pseudoalignments for the reads ... done
## [quant] processed 24,729,208 reads, 20,253,657 reads pseudoaligned
## [quant] estimated average fragment length: 174.106
## [   em] quantifying the abundances ... done
## [   em] the Expectation-Maximization algorithm ran for 1,069 rounds
## [bstrp] running EM for the bootstrap: 100
## ...

kallisto はマッピングを伴わないので、BAM/SAM ファイルは生成されない。その代わり、推定された発現量はテキストファイル abundances.tsv に保存される。また、-b によるブートストラップの結果は abundances.h5 ファイルに保存される。abundances.h5 は HDF5 フォーマットで、sleuth に発現変動遺伝子検出時の入力するファイルである。また、kallisto を実行した時のオプションなどは run_info.json に保存される。

ls SRR5383650_exp_kallisto
## abundance.h5  abundance.tsv  run_info.json

head SRR5383650_exp_kallisto/abundance.tsv
## target_id	length	eff_length	est_counts	tpm
## TraesCS1A02G000100.1	1308	1134.89	13.5027	0.619261
## TraesCS1A02G000200.1	1416	1242.89	481.736	20.1736
## TraesCS1A02G000300.1	795	621.894	39.1928	3.28019
## TraesCS1A02G000400.1	3055	2881.89	1204.51	21.7541
## TraesCS1A02G000500.1	408	235.18	1.02046	0.225842
## TraesCS1A02G000600.1	297	126.661	0	0
## TraesCS1A02G000700.1	270	102.135	0	0
## TraesCS1A02G000800.1	378	205.287	0	0
## TraesCS1A02G000900.1	2987	2813.89	1147.03	21.2167

run_info.json を確認すると、81.9% が擬似アラインメントされて、56.2% が uniquely であった。pseudoaligned の割合が高いにもかかわらず、uniquely pseudoaligned の割合が非常に低い。これは小麦が 6 倍体であり同一遺伝子のコピー数が非常に多いのが一つの原因となっているかもれない。kallisto のインデックスを作成するとき、k-mer のサイズを長くすることで、uniquely pseudoaligned を減らせるかもしれない。ただ、kallisto では最長で 31 bp しか指定できないようである。

cat SRR5383650_exp_kallisto/run_info.json
## {
## 	"n_targets": 133744,
## 	"n_bootstraps": 100,
## 	"n_processed": 24729208,
## 	"n_pseudoaligned": 20253657,
## 	"n_unique": 13907889,
## 	"p_pseudoaligned": 81.9,
## 	"p_unique": 56.2,
## 	"kallisto_version": "0.44.0",
## 	"index_version": 10,
## 	"start_time": "Mon Aug 27 11:39:54 2018",
## 	"call": "kallisto quant -i IGWSCv1_kallisto -o SRR5383650_exp_kallisto -b 100 SRR5383650_1.clean.fastq.gz ## SRR5383650_2.clean.fastq.gz"
## }

sleuth で発現変動転写産物（遺伝子）を検出したい場合は、R の sleuth パッケージのマニュアルを参照して、kallisto の結果を入力として与えればいい。一方で、edgeR や DESEq2 などを使いたい場合は、R の tximport パッケージを介して kallisto の結果（abundances.tsv または abundances.h5）を edgeR や DESeq2 に入力するのが便利（実行例）。

References

Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016, 34(5):525-7. DOI: 10.1038/nbt.3519
IWGSC. Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science. 2016, 361(6403):eaar7191. DOI: 10.1126/science.aar7191
Ma J, Li R, Wang H, Li D, Wang X, Zhang Y, Zhen W, Duan H, Yan G, Li Y. Transcriptomics Analyses Reveal Wheat Responses to Drought Stress during Reproductive Stages under Field Conditions. Front Plant Sci. 2017, 8:592. DOI: 10.3389/fpls.2017.00592