macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

真核生物のRNA-seqデータ解析

今回は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

f:id:kazumaxneo:20170324214520p:plain

となっている。

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を使うワークフロー。