今回はhumanのデータを使う。重いのでchromosome19に限定して解析する。
------------準備------------
1、データ
fastaファイル。 FTP Downloadからhumanのfastaファイルをダウンロードする。
Homo_sapiens.GRCh38.dna_rm.chromosome.19.fa
rRNAのfastaファイル。ここではgitに上がっているのをダウンロード。
IlluminaPE/humRibosomal.fa at master · Magdoll/IlluminaPE · GitHub
配列をコピペしてrRNA.faとして保存。
gtfファイル。 FTP Downloadからhumanのgtfファイルをダウンロードする。
Homo_sapiens.GRCh38.87.gtf
fastqファイル 。http://trace.ddbj.nig.ac.jp/DRASearch/experiment?acc=DRX015048に行って、右端のfastqを選択。ウィンドウを開くか聞いてくるのでゲストを選択して開く。
DRR016695_1.fastq.bz2とDRR016695_2.fastq.bz2を解析フォルダにコピペしてダウンロード(bz2圧縮で2GBx2)。
2、解析プログラム
homebrewで実行環境を整える。homebrewはこちら
brew install tophat
brew install bowtie2
brew install samtools
brew install cufflinks
------------解析------------
はじめにrRNAにmapされたリードを除く。それからgtfに記載されている既存exsonにマッピングして、fpkmを算出する。
1、rRNA.faのindex作成 (bowit2を使う)
bowtie2-buid humRibosomal.fa humRibosomal
2、chr19.faのindex作成 (bowit2を使う)
bowtie2-buid Homo_sapiens.GRCh38.dna_rm.chromosome.19.fa chr19
3、rRNAへのマッピング解析 (tophatを使う)
bowtie2 -p 8 --un DRR016695_1_without_rRNA.fastq -x humRibosomal -U DRR016695_1.fastq -S DRR016695_1_rRNA_mapped.sam
-p スレッド数
--un <path> write unpaired reads that didn't align to <path>
-x rRNAインデックス作成時の名前
-U インプットに使うunpairのfastq
-S マッピング結果出力名
もう片方のfastqも同様に処理する。
bowtie2 -p 8 --un DRR016695_2_without_rRNA.fastq -x humRibosomal -U DRR016695_2.fastq -S DRR016695_2_rRNA_mapped.sam
rRNAにマッピングされなかったリードとして出力されたDRR016695_1_without_rRNA.fastqとDRR016695_2_without_rRNA.fastqを以降の解析に使う。
4、chr19へのマッピング解析
tophat -G Homo_sapiens.GRCh38.87.gtf -p 8 -o human_output --no-novel-juncs chr19 DRR016695_1_without_rRNA.fastq DRR016695_2_without_rRNA.fastq
-p たぶんスレッド数。増やしすぎるとエラーになる。CPUが多くても10以内に止める。
-o 出力フォルダ名
-G エキソンの識別に使うgtfファイル。上でダウンロードしたファイルを使う。
-no-novel-juncs アノテーションにない新規のジャンクションを破棄
chr19 上のbowtie2のステップで作成したインデックス名。
-Gの代わりに-gを使うと、既存のエキソン情報を使いながら、新規スプライシング産物も探索するようになる。
解析が終わると指定した出力フォルダの中にaccepted.bamができる。これを次の定量に使う。
5、fpkmの算出
cufflinks -G Homo_sapiens.GRCh38.87.gtf -N --compatible-hits-norm -M rRNA.gtf -o output_fpkm human_output/accepted_hits.bam
-o 出力フォルダ名
-G エキソンの識別に使うgtfファイル。
-N --compatible-hits-norm 既存のアイソフォームにマッピングされるものだけでFPKMを計算
-M rRNA.gtfにマッピングされるリードを除く。
解析が終わると
genes.fpkm_tracking
isoforms.fpkm_tracking
transcripts.gtf
が指定フォルダに出力される。gene.trackは遺伝子単位の発現量、isoform.trackは愛想フォーム単位の発現量のファイル。中身を見ると
>head -3 pkm human/ genes.fpkm_tracking
となっている。
tracking id 遺伝子ID。
class code アイソフォームのタイプ。基本的に-
nearest ref id 近傍の遺伝子ID。
gene id 遺伝子ID(1列目と同じ)。
gene short name 遺伝子名。
tss id 転写開始点ID。
locus ポジション。
length 遺伝子の全長。
coverage カバレッジ。
FPKM サンプルFPKM。
FPKM conf lo 下限95%区間のサンプルFPKM。
FPKM conf hi 上限95%区間のサンプルFPKM。
FPKM status サンプルのFPKMステータス。OKかLOW、Hi、FAIL。
複数サンプルがある場合は
q0_FPK(サンプル1) q1_FPKM(サンプル2)などの名前で区別される。
一般的な 全体のフローは以下で説明しています。
featurecountsとedgeRを使うワークフロー。