ペアエンド の RNA-Seq データ解析のパイプラインをメモします。
はじめに
Linux OS 上で Shellscript を用いて各種NGS解析ツールを走らせて解析をします。解析内容は以下です。
- md5sumによるデータの確認
- Trimmomaticによるfastqデータのクオリティチェック
- HISAT2とsamtoolsによるマッピング
- +/-鎖にbamファイルを分割
- bamファイルのソートとインデックス作成
- depthファイルの作成
- featureCountsによるカウントデータの作成
データの確認
NGSによるシーケンスを外注した場合、fastqデータと一緒にmd5sumの値が納品されます。
自身で納品されたfastqデータのmd5sumを計算し、それが納品されたmd5sumの値と等しいか確認します。
もし等しければ、納品されたfastqデータに破損等がないと判断します。
まずは以下のようにして得られたfastqのmd5sumを計算します。
$ pwd /home/neko/data $ md5sum `ls *.fastq.gz` > md5sum.txt
md5sum.txt は以下のような内容になります。
2df0b1f2766fcade5484ab21b80b7306 sample1.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample2.fastq.gz 0cfe18d798693c93abc86f2f8a675483 sample3.fastq.gz 25a6fab8c3cf283d2501c3f1f6bd2811 sample4.fastq.gz 8d06752557e09a7e7e5ffb5dba96e225 sample5.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample6.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample7.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample8.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample9.fastq.gz
次にdiffコマンドを用いてデータが正しいか調べます。-yオプションで左右並べて表示できます。
$ diff -y md5sum.txt md5sum_true.txt > diff.txt
md5sum.txt と md5sum_true.txt の各行が同じである場合、上記の diff コマンドの出力(diff.txt)は以下のようになります。
2df0b1f2766fcade5484ab21b80b7306 sample1.fastq.gz 2df0b1f2766fcade5484ab21b80b7306 sample1.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample2.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample2.fastq.gz 0cfe18d798693c93abc86f2f8a675483 sample3.fastq.gz 0cfe18d798693c93abc86f2f8a675483 sample3.fastq.gz 25a6fab8c3cf283d2501c3f1f6bd2811 sample4.fastq.gz 25a6fab8c3cf283d2501c3f1f6bd2811 sample4.fastq.gz 8d06752557e09a7e7e5ffb5dba96e225 sample5.fastq.gz 8d06752557e09a7e7e5ffb5dba96e225 sample5.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample6.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample6.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample7.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample7.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample8.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample8.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample9.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample9.fastq.gz
もし不一致な行があった場合、例えば以下のように出力されます。| や > などの記号が入ります。以下では下から2行目に不一致な箇所があることが分かります。
2df0b1f2766fcade5484ab21b80b7306 sample1.fastq.gz 2df0b1f2766fcade5484ab21b80b7306 sample1.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample2.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample2.fastq.gz 0cfe18d798693c93abc86f2f8a675483 sample3.fastq.gz 0cfe18d798693c93abc86f2f8a675483 sample3.fastq.gz 25a6fab8c3cf283d2501c3f1f6bd2811 sample4.fastq.gz 25a6fab8c3cf283d2501c3f1f6bd2811 sample4.fastq.gz 8d06752557e09a7e7e5ffb5dba96e225 sample5.fastq.gz 8d06752557e09a7e7e5ffb5dba96e225 sample5.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample6.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample6.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample7.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample7.fastq.gz dab3d4103742aa1e4acd92c87f8f70c2 sample8.fastq.gz | dab3d4103742aa1e4acd92c87f8f70c1 sample8.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample8.fastq.gz 69a42ad1e16d1fcb71753ba1016c345d sample9.fastq.gz
上記のようにして納品されたデータに間違いがないことを確認できたら、以下のようにして解析フォルダへfastqファイルを移動させ、次のステップに進みます。
$ mkdir -p /home/neko/rnaseq $ mv *.fastq.gz /home/neko/rnaseq
クオリティチェック
Trimmomatic は fastq データのクオリティを調べて、クオリティが低いリードや NGS のアダプター配列を除去します。
以下のようなシェルスクリプト(trimmo.sh)を/home/neko/rnaseq/scriptに用意します。
BASE=/home/neko/rnaseq
fname=$1
cd $BASE
Cpu=14
# 1. Quality Filtering (Trimmomatic)
Min=36
Read=$BASE/fastq
Dir=$BASE/fastq_clean
Rep=$BASE/reports/main/trimmo
Trimmo=/lustre/home/70834888/bin/Trimmomatic-0.39/trimmomatic-0.39.jar
Iclip=/lustre/home/70834888/bin/Trimmomatic-0.39/adapters/all.SE.fa
mkdir -p $Dir $Rep
java -jar $Trimmo \
PE -phred33 -threads $Cpu \
-trimlog $Rep/$fname.log \
$Read/${fname}_1.fq.gz \
$Read/${fname}_2.fq.gz \
$Dir/${fname}_paired_1.fq \
$Dir/${fname}_unpaired_1.fq \
$Dir/${fname}_paired_2.fq \
$Dir/${fname}_unpaired_2.fq \
ILLUMINACLIP:$Iclip:2:30:10 \
MINLEN:$Min上記のスクリプトで、$BASE と $Trimmo と $Iclip に代入するパスはそれぞれの環境に合わせて変えて下さい。
なお、Trimmomatic のインストールはこちらを参考にして下さい。
上記のスクリプトが用意できたら、以下のようにして実行して下さい。
$ bash trimmo.sh sample1.fastq
マッピング
本記事ではマッピングのツールに HISAT2 を用います。HISAT2 のインストールはこちら。
マッピングにはレファレンス配列(ゲノム配列など)のファイルと、それに対するインデックスファイルの作成が必要になります。
レファレンスファイルを用意するには様々な方法がありますが、例えば Ensembl 等からダウンロードできます。ヒトのレファレンスファイルとなるゲノムのファイルはこちらの”Homo_sapiens.GRCh38.dna.toplevel.fa.gz”というファイルをダウンロードすればOKです。解凍は以下のようにして行います。
$ gzip -d Homo_sapiens.GRCh38.dna.toplevel.fa.gz
インデックスファイルの作成は、HISAT2 を用いて以下のようにできます。環境にも依りますが、30分以上の時間がかかります。
Hisat2_build=/home/neko/bin/hisat2-2.2.1/hisat2-build Ref=/home/neko/src/Homo_sapiens.GRCh38.dna.toplevel.fa Index_name=/home/neko/src/Homo_sapiens.GRCh38.dna.toplevel $Hisat2_build $Ref $Index_name
インデックスの作成が完了したら、以下のようにマッピングのスクリプト(mapping.sh)を用意します。
# 2. Mapping (HISAT2, samtools)
Read=$BASE/fastq_clean
Dir=$BASE/bam
Hisat2=/home/neko/bin/hisat2-2.2.1/hisat2
Ref=/home/neko/src/Homo_sapiens.GRCh38.dna.toplevel.fa
IndexName=/home/neko/src/Homo_sapiens.GRCh38.dna.toplevel
samtools=/lustre/home/70834888/bin/samtools/samtools
mkdir -p $Dir
$Hisat2 -x $IndexName \
-1 $Read/${fname}_paired_1.fq \
-2 $Read/${fname}_paired_2.fq \
-p $Cpu | $samtools view -Shb - > $Dir/$fname.bam$BASE と $Cpu と $fname は trimmo.sh と同じです。
samtoolsのインストール方法はこちら。
上記のスクリプトが用意できたら、以下のようにして実行します。
$ bash mapping.sh
+/-鎖へのbamファイルの分割
RNA-Seq のライブラリ調製方法によっては、リードがゲノムの+鎖と-鎖のどちらにマッピングされたか識別可能です。
RNA-Seqの解析では必ずしも実施する必要はないですが、以下のようにして+鎖と-鎖への bamファイルの分割ができます。
$samtools view -b -f 0 $Dir/$fname.bam > $Dir/$fname.plus.bam $samtools view -b -f 16 $Dir/$fname.bam > $Dir/$fname.minus.bam
$samtools、$Dir、$fname は mapping.sh と同じです。
bamファイルのソートとインデックス作成
マッピング後の多くのプロセスでは、bamファイルがソートされていることが必要です。また、bamファイルのインデックスファイルが必要なことも多いです。
そこで、マッピングした後は、通常bamファイルのソートとインデックス作成を以下のようにして行います。
$samtools sort -@ $Cpu -o $Dir/$fname.plus.refsort.bam $Dir/$fname.plus.bam $samtools index $Dir/$fname.plus.refsort.bam $samtools sort -@ $Cpu -o $Dir/$fname.minus.refsort.bam $Dir/$fname.minus.bam $samtools index $Dir/$fname.minus.refsort.bam
$samtools、$Cpu、$Dir、$fname は mapping.sh と同じです。
depthファイルの作成
bamファイルをソートしてインデックスを作成した後の解析は、目的によって様々に分かれていきます。
この節では depth ファイルを作成する方法を示します。
depth ファイルはゲノム上の各位置に対してリードが何本マッピングされたかを示すファイルです。
まずは、以下のような、depth.sh を作成します。
# 3. Depth (samtools) Read=$BASE/bam Dir=$BASE/depth mkdir -p $Dir $samtools depth $Read/$fname.plus.refsort.bam > $Dir/$fname.plus.depth $samtools depth $Read/$fname.minus.refsort.bam > $Dir/$fname.minus.depth
$BASE、$samtools、$fname は mapping.sh と同じです。
上記のスクリプトを以下のように実行すれば、depth ファイルはできます。
$ bash depth.sh
カウントデータの作成
この節ではソートされた bamファイルを入力に、featureCountsを用いてカウントデータを作成する方法を示します。
featureCounts は、マッピング結果からどの feature 領域（gene、CDS、exon など）に何本のリードがマッピングされたのか集計するプログラムで、ゲノム上のどの領域に feature があるか記載されたGTFファイルが必要です。
GTFファイルはマッピングの際にレファレンスファイルをダウンロードした Endembl からダウンロードできます。ヒトではこちらの”Homo_sapiens.GRCh38.111.gtf.gz”というファイルをダウンロードすればOKです。解凍は以下のように行います。
$ gzip -d Homo_sapiens.GRCh38.111.gtf.gz
カウントについては、以下のように、エクソン領域にマッピングされたリードのみを計上し（-t exon）、その後、遺伝子名ごとに集計します（-g gene_id）。エクソン領域や遺伝子名のアノテーションファイルは -a の引数として与えます。
# 4. Counts(featureCounts)
Read=$BASE/bam
Dir=$BASE/count
featureCounts=/home/neko/bin/featureCounts
mkdir -p $Dir
$featureCounts -t exon \
-g gene_id \
-a annotation.gtf \
-o counts.txt mapping_results.bam$BASEは mapping.sh と同じです。
コメント