ART: シミュレーション用のRNA-seqデータを作成する

NGS RNA-seq Simulation ART

テスト用のFASTQファイルの作成

ART (a next-generation sequencing read simulator)を用いて、テスト用のFASTQファイルを作成します。

ARTのインストール

※ Biocondaをインストールしていることを前提にしています。
ARTをインストールします。

$ conda intall art

ARTが正しくインストールされたか確認します。

$ art_illumina --help

RefseqのGTFファイルをダウンロード

ARTを使ってテスト用のFASTQファイルを作るのですが、今回はRNA-seqのデータを作ろうと思います。

そこで、元となる配列を用意したいと思います。 TranscriptomeのデータをGTFファイルと取得し、それをFASTQのファイルに変換したいと思います。

GencodeのGTFファイルでもいいですが、ここでは必要なTranscriptomeデータが含まれておりコンパクトなRefseqのデータを使います。

illuminaのiGenomeサイトから、hg38のデータを丸ごとダウンロードします。この中に、RefSeqのGTFファイルも含まれています。

image.png (182.7 kB)

ダウンロードしたら、解凍します。

$ tar zxvf Home_sapiens_UCSC_hg38.tar.gz

解凍が終わったら、Homo_sapiensというディレクトリができるはずです。 Homo_sapiens/UCSC/hg38/Annotation/Archives/archive-2015-08-14-08-18-15/Genes/にあるgenes.gtfというファイルがRefSeqのGTFファイルになります（ちょっと古いデータですが、実際の解析に用いるわけではないので、良しとしましょう）。

これだけ適当なディレクトリに移動もしくはコピーしておき、それ以外のファイルは削除します。

このままだと、何のファイルかわからないのでリネームしておきます。

$ mv genes.gtf RefSeq_hg38_2015-08-19.gtf

GTFファイルからBEDファイルを作成

以下のスクリプトを使って、GTFファイルをBEDファイルに変換する。 gtf2bed.pl (3.8 kB)

$ perl gtf2bed.pl RefSeq_hg38_2015-08-19.gtf > RefSeq_hg38_2015-08-19.bed

このままだと、データが大きすぎるのでmRNAに限定し、かつRepresentative transcriptsをできるだけ抽出するようにします。
ExtractMrna.py (360 B)

$ python ExtractMrna.py RefSeq_hg38_2015-08-19.bed RefSeq_hg38_2015-08-19_NM_slim.bed

ヒトゲノム配列をダウンロード

image.png (1.5 MB)

image.png (1.5 MB)

image.png (300.8 kB)

image.png (71.1 kB)

image.png (273.5 kB)

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
$ gunzip hg38.fa.gz

BEDファイルからFASTAファイルを作成

BEDファイルからFASTA ファイルへ変換するために、Bedtoolsを使います。 まず、Bioconda経由でBedtoolsをインストールします。

$ conda install bedtools

BEDファイルをFASTAファイルに変換します。 ここで、先ほどダウンロードしたhg38.faファイルを使います。

$ bedtools getfasta -name -s -split -fi hg38.fa -bed RefSeq_hg38_2015-08-19_NM_slim.bed > RefSeq_hg38_2015-08-19_NM_slim.fa

-name: BEDファイルのname列の値を各配列の名称とする。
-s: ゲノム上での向き (Strandness)を考慮する（Transcriptomeなのでもちろん考慮します）。
-split: 12BED (12列のデータがあるBEDふぁいる)がInputである場合に指定 。
-fi: Referenceになる配列（ゲノム配列）のFASTAファイルを指定。
-bed: TranscriptomeのBEDファイルを指定。

このコマンドの実行には、FASTAファイルのインデックスファイルが必要になりますが、ファイルがない場合、自動的に生成してくれるので問題ないです。

> index file hg38.fa.fai not found, generating...

生成されたTranscriptomeのFATSAファイルは、標準出力で出力されます。

TranscriptomeのFASTAファイルからFASTQファイルを生成

今回は、illuminaのFASTQシーケンスデータ (RNA-seq, single-end)を生成したいと思います。ここでは、HiSeq2000のデータを作ります。

$ art_illumina -ss HS20 -i RefSeq_hg38_2015-08-19_NM_slim.fa -o RefSeq_hg38_2015-08-19_NM_slim -l 36 -f 2 --noALN
$ art_illumina -ss HS20 -i RefSeq_hg38_2015-08-19_NM_slim.fa -o RefSeq_hg38_2015-08-19_NM_slim_paired -l 36 -f 1 -p -m 300 -s 10 --noALN

-i --in: fasta形式の入力ファイルの指定。
-o --out: fastq形式の出力ファイルの指定 (拡張子はつけなくて良い)。
-ss --seqSys: シーケンサーの種類を選択。
指定する引数(System ID) - シーケンサー名 (リード長)

  • GA1 - GenomeAnalyzer I (36bp,44bp),
  • GA2 - GenomeAnalyzer II (50bp, 75bp)
  • HS10 - HiSeq 1000 (100bp)
  • HS20 - HiSeq 2000 (100bp)
  • HS25 - HiSeq 2500 (125bp, 150bp)
  • HS10 - HiSeq 1000 (100bp)
  • HS20 - HiSeq 2000 (100bp)
  • HS25 - HiSeq 2500 (125bp, 150bp)
  • HSXn - HiSeqX PCR free (150bp)
  • HSXt - HiSeqX TruSeq (150bp)
  • MinS - MiniSeq TruSeq (50bp)
  • MSv1 - MiSeq v1 (250bp)
  • MSv3 - MiSeq v3 (250bp)
  • NS50 - NextSeq500 v2 (75bp)

-l --len: リード長の指定。
-f --fcov: リードカバレージの指定。
-p --paired: Paired-endのデータとして出力（指定しない場合、Single-endのデータとして出力）。
m --mflen: フラグメントの平均サイズの指定（Paired-endのときのみ）。
-s --sdev: フラグメントのサイズの標準偏差の指定（Paired-endのときのみ）。
-sam --samout: SAMファイルを出力。
-na --noALN: アライメントファイル (ALN)を出力しない。
別にリード長は、指定された長さよりも短くても大丈夫です。

もっとファイルサイズを小さくしたい場合は、以下のように特定の行を抽出したファイルを作成します。

$ sed -n 1,200000p RefSeq_hg38_2015-08-19_NM_slim.fq > RefSeq_hg38_2015-08-19_NM_ultraSlim.fq