#※記事を読んでくださるみなさまへ
・対象はプログラミングもバイオインフォ解析も全く未経験のバイオ研究者です
・筆者は基礎研究を初めて約1年3ヶ月の大学院生です
・いわゆる”DRY”な研究を始めたのは約5ヶ月前からです
・そのため記事の内容にはとんでもない間違いやずっと良い別の方法が存在するかもしれません
・また、筆者と同じレベルの人が0から始めるときの参考になるよう、説明は結構まだるっこいです
・ご指摘の度に記事をブラッシュアップしていきたいです、よろしくお願いします
過去記事一覧
①環境構築まで
②NGSデータのダウンロードとクオリティチェック
③NGSデータのトリミング
④Pipeline (hisat2, samtools, stringtie)
⑤番外編：コードのブラッシュアップ
⑥Rによる発現解析
前回の続きから：残りのfastq.gzファイルもトリミング
BoWのタイトルバーを右クリック、プロパティから簡易編集モードをOnにしてください。lsで表示したファイル名をコピペして、メモ帳などに貼っておきましょう。
そしたらミスをしないように、Wordに前回のコードを貼り付けて、数字の部分を検索・置換して各サンプルデータ用のコードを作りましょう（もっとうまいやり方がわかる方いたら教えてください）。
#HISAT2でのマッピング
見づらいので使うファイルだけディレクトリを作って移してしまいましょう（オプション）。ついでに使わないunpairedのファイルは消してしまいます。
mkdir paired
mv *_paired.fq.gz paired
rm *_unpaired.fq.gz
# ＊はワイルドカード。この場合末尾だけが合致するもの全て、ということ
以下の解析ではCPUの「スレッド数」を知ることが大事になります。要するに使用できるCPUのパワー上限みたいなものです。
lscpu
Architecture: x86_64
CPU 操作モード: 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
コアあたりのスレッド数:2
ソケットあたりのコア数:2
Socket(s): 1
ベンダー ID: GenuineIntel
CPU ファミリー: 6
モデル: 78
ステッピング: 3
CPU MHz: 2208.000
BogoMIPS: 4416.00
仮想化: VT-x
みたいな表示が出て、4✕2＝8スレッド使用可能だとわかります。いよいよマッピングしましょう。
cd chrX_data/
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188044_chrX_1_paired.fq.gz -2 samples/paired/ERR188044_chrX_2_paired.fq.gz -S ERR188044_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188104_chrX_1_paired.fq.gz -2 samples/paired/ERR188104_chrX_2_paired.fq.gz -S ERR188104_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188234_chrX_1_paired.fq.gz -2 samples/paired/ERR188234_chrX_2_paired.fq.gz -S ERR188234_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188245_chrX_1_paired.fq.gz -2 samples/paired/ERR188245_chrX_2_paired.fq.gz -S ERR188245_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188257_chrX_1_paired.fq.gz -2 samples/paired/ERR188257_chrX_2_paired.fq.gz -S ERR188257_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188273_chrX_1_paired.fq.gz -2 samples/paired/ERR188273_chrX_2_paired.fq.gz -S ERR188273_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188337_chrX_1_paired.fq.gz -2 samples/paired/ERR188337_chrX_2_paired.fq.gz -S ERR188337_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188383_chrX_1_paired.fq.gz -2 samples/paired/ERR188383_chrX_2_paired.fq.gz -S ERR188383_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188401_chrX_1_paired.fq.gz -2 samples/paired/ERR188401_chrX_2_paired.fq.gz -S ERR188401_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188428_chrX_1_paired.fq.gz -2 samples/paired/ERR188428_chrX_2_paired.fq.gz -S ERR188428_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR188454_chrX_1_paired.fq.gz -2 samples/paired/ERR188454_chrX_2_paired.fq.gz -S ERR188454_chrX.sam
hisat2 -p 8 --dta -x indexes/chrX_tran -1 samples/paired/ERR204916_chrX_1_paired.fq.gz -2 samples/paired/ERR204916_chrX_2_paired.fq.gz -S ERR204916_chrX.sam
各オプションについて解説すると、-pは使用するスレッド数（デフォルトは1）、--dtaは以下の解析で使用するのがStringtieのときのオプション、-xは「インデックスファイル」の場所。インデックスファイルは、高速化のための目次ファイルのことらしいです。主な種のIndexファイルはHISATのWebページにあります。-1、-2はForward, Revereseリードファイル（多分順番はどっちでもよさそう）。-Sは出力されるファイルのファイル名です。
どうもマッピングエラ－がなくなるせいか、トリミングしてからのほうがしないでやるよりマッピング速い気がします。
（5/8追記）
シングルエンドリードのときのコードを書き忘れていました。
hisat2 -p 8 --dta -x （インデックスファイルのあるディレクトリパス）/（インデックスファイルの.数字.ht2の前の部分の文字列） -U （インプットファイルのディレクトリパス）/ファイル名 -S （保存名）.sam
-Uがシングルエンドのときのインプットファイルのオプション。インデックスファイルを替えたとき（マウス→ヒトとか、ゲノムリファレンス→トランスクリプトームリファレンスみたいなとき）にハマったのだが、-xオプションはペアエンドであってもシングルエンドであっても、上記のような特殊な記載法になる。パスで終わりになるのではなく、インデックスファイルの頭文字になる点が要注意。
#SAMファイルをBAMファイルに変換
以下の解析を高速化するため、SAMファイルという人間がわかる形式のファイルを、BAMファイルという、機械が理解しやすい「バイナリファイル」に変換します。
samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
samtools sort -@ 8 -o ERR188104_chrX.bam ERR188104_chrX.sam
samtools sort -@ 8 -o ERR188234_chrX.bam ERR188234_chrX.sam
samtools sort -@ 8 -o ERR188245_chrX.bam ERR188245_chrX.sam
samtools sort -@ 8 -o ERR188257_chrX.bam ERR188257_chrX.sam
samtools sort -@ 8 -o ERR188273_chrX.bam ERR188273_chrX.sam
samtools sort -@ 8 -o ERR188337_chrX.bam ERR188337_chrX.sam
samtools sort -@ 8 -o ERR188383_chrX.bam ERR188383_chrX.sam
samtools sort -@ 8 -o ERR188401_chrX.bam ERR188401_chrX.sam
samtools sort -@ 8 -o ERR188428_chrX.bam ERR188428_chrX.sam
samtools sort -@ 8 -o ERR188454_chrX.bam ERR188454_chrX.sam
samtools sort -@ 8 -o ERR204916_chrX.bam ERR204916_chrX.sam
sortはsamtoolsの機能のうち、並べ替え機能。何で並べ替えているのかはWebページに書いてあります。”leftmost”と書いてあるので、多分転写産物名です。-@はスレッドの使用数（デフォルトは1）。-oがBAM形式で出力するよ、という指示で、拡張子だけ変えた名前にしておくと良いでしょう。一番お尻が入力ファイルです。
#サンプルごとのGTFファイルの作成
ではこれらのBAMファイルからGTFファイルという、遺伝子アノテーション情報を付加したファイルを作成しましょう。
stringtie -p 8 -G genes/chrX.gtf -o ERR188044_chrX.gtf -l ERR188044 ERR188044_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188104_chrX.gtf -l ERR188104 ERR188104_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188234_chrX.gtf -l ERR188234 ERR188234_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188245_chrX.gtf -l ERR188245 ERR188245_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188257_chrX.gtf -l ERR188257 ERR188257_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188273_chrX.gtf -l ERR188273 ERR188273_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188337_chrX.gtf -l ERR188337 ERR188337_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188383_chrX.gtf -l ERR188383 ERR188383_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188401_chrX.gtf -l ERR188401 ERR188401_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188428_chrX.gtf -l ERR188428 ERR188428_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR188454_chrX.gtf -l ERR188454 ERR188454_chrX.bam
stringtie -p 8 -G genes/chrX.gtf -o ERR204916_chrX.gtf -l ERR204916 ERR204916_chrX.bam
-pは使用するスレッド数（デフォルト1）。-Gはリファレンスアノテーションファイル（遺伝子アノテーション情報が載っている）。-oは出力するファイル名。-lはラベル名で、出力ファイルの転写産物名の頭につけられるラベル。たぶんこれをサンプルの区別に使っている。お尻が入力するファイル。
（5/9追記）
各動物種のアノテーションファイルはここから落とせるみたいです。マウスだとこれのMus_musculus.GRCm38.88.chr_patch_hapl_scaffであってるんでしょうか？わかる方いたら教えてください。
#GTFファイルの結合
Stringtieというプログラムは、サンプルごとに得られたトランスクリプト情報とそれをマージしたトランスクリプト情報を突き合わせることで発現量の再推定を行い、正確性を上げているそうです。以降はそれに必要なマージファイルを作ります。コードは簡単です。
stringtie --merge -p 8 -G genes/chrX.gtf -o stringtie_merged.gtf mergelist.txt
--mergeはマージするよ、というオプション。-pは使用するスレッド数（デフォルトは1）。-Gはリファレンスとなる遺伝子アノテーション情報が記載されたファイル。-oは出力するファイル名。mergelist.txtは自前のデータでやるときに一番ミスりやすいのですが、入力ファイル一覧が記載されたtxtファイルです。cat mergelist.txtとすると以下のような内容になっています。
ERR188044_chrX.gtf
ERR188104_chrX.gtf
ERR188234_chrX.gtf
ERR188245_chrX.gtf
ERR188257_chrX.gtf
ERR188273_chrX.gtf
ERR188337_chrX.gtf
ERR188383_chrX.gtf
ERR188401_chrX.gtf
ERR188428_chrX.gtf
ERR188454_chrX.gtf
ERR204916_chrX.gtf
なので自前のデータでの解析のときには自力で作りましょう。
#Ballgown用データセットの作成
いよいよ最後です。実際にRのballgownパッケージで解析するためのデータセットを作ります。
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188044/ERR188044_chrX.gtf ERR188044_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188104/ERR188104_chrX.gtf ERR188104_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188234/ERR188234_chrX.gtf ERR188234_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188245/ERR188245_chrX.gtf ERR188245_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188257/ERR188257_chrX.gtf ERR188257_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188273/ERR188273_chrX.gtf ERR188273_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188337/ERR188337_chrX.gtf ERR188337_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188383/ERR188383_chrX.gtf ERR188383_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188401/ERR188401_chrX.gtf ERR188401_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188428/ERR188428_chrX.gtf ERR188428_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188454/ERR188454_chrX.gtf ERR188454_chrX.bam
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR204916/ERR204916_chrX.gtf ERR204916_chrX.bam
-eは-Gで与えたリファレンスアノテーション以外の除去。-Bは以下の解析で使うballgown用のファイルを出力するような指示。-pは使用するスレッド数（デフォルトは1）。-oは出⼒するファイル名。後のballgown()関数の引数でサンプルごとのディレクトリを指定するので、
ballgownというディレクトリを作成し、サンプルごとのディレクトリを作り、その中に入れる形になっています。また、実際にはgtfファイルだけではなく、先に指示したballgown用の.ctabといったデータも出力されるみたいです。お尻は ⼊⼒するファイル。
長かったですが、以上でRでの解析までの下準備が全て終わりました。次回以降はRStudioを使って、サンプルごとの発現解析を行っていこうと思います。