2017-03-24

Artemis使用法

2020 3/1 リンク追加

Artemisはサンガー研で開発されたゲノムデータをビジュアルで見れるソフト。ゲノムのデータ（genebakファイルなど）を読み込ませると、遺伝子単位まで拡大して確認できる。ダウンロードはここから http://www.sanger.ac.uk/science/tools/artemis

f:id:kazumaxneo:20170324171301j:plain

Macなら、写真のMac OSをクリック。自動的にダウンロードが開始される。ダウンロードした.dmgファイルを解凍するとArtemisが表示される。Applicationにコピーしてインストール完了。起動したらFile -> open... を選択し、

f:id:kazumaxneo:20170324171342j:plain

genebakフォーマットのファイル(.gbk)を選択する。

読み込むと以下のようにorfを可視化して表示できる。

f:id:kazumaxneo:20170324171418j:plain

拡大縮小方法はやや特殊。

上の右端に小さな灰色ゲージが見えるが、これをマウスで上下に引っ張れば拡大縮小される。ゲノム全体を表示することもできる。

f:id:kazumaxneo:20170324171451j:plain

まとめ直しました。

2017-03-24

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

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）。

２、解析プログラム

homebrewで実行環境を整える。homebrewはこちら

brew install tophat

brew install bowtie2

brew install samtools

brew install cufflinks

------------解析------------

はじめにrRNAにmapされたリードを除く。それからgtfに記載されている既存exsonにマッピングして、fpkmを算出する。

１、rRNA.faのindex作成 (bowit2を使う）

bowtie2-buid humRibosomal.fa humRibosomal

２、chr19.faのindex作成 (bowit2を使う）

bowtie2-buid Homo_sapiens.GRCh38.dna_rm.chromosome.19.fa chr19

３、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を以降の解析に使う。

４、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（１列目と同じ）。

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（サンプル１）　q1_FPKM（サンプル２）などの名前で区別される。

一般的な全体のフローは以下で説明しています。

featurecountsとedgeRを使うワークフロー。

2017-03-24

ggplot2によるグラフ作成

other

Contig_Hunter_version_0.1.plで解析すると、-R_ggplot2というフォルダの中に４つのファイルが出力される。それをRに読み込ませる。*1

coverage.txt GC.txt

length.txt name.txt

Rのワーキングディレクトリに４つのファイルをコピーし、Rのターミナル環境で以下のコマンドを実行（ggplot2とscalesライブラリがインストールされてない場合は、予めインストールしておくこと）。

library(ggplot2)

library(scales)

( x <- read.table("coverage.txt", header=T) ) #カバレッジ

( y <- read.table("length.txt", header=T) ) #x軸も同様に入力

( z <- read.table("GC.txt", header=T) ) #3番目の要素

df <- data.frame(x=x, y=y,z)

g <- ggplot(df,aes(x = x,y = y,z))

g <- g + xlab("coverage") # x 軸ラベル

g <- g + ylab("length") # y 軸ラベル

#グラフ6 目盛りを統一

xbreaks <- c(10,100,500)

ybreaks <- seq(100,10000,1000000)

#最後に以下の１行をペースト

g + geom_point(aes(colour=z,size = GC),alpha = 0.5) + scale_size_area(name = "GC", max_size = 5) + scale_x_log10(labels=trans_format("log10",math_format(10^.x))) + scale_y_log10(labels=trans_format("log10",math_format(10^.x))) + scale_colour_gradientn(colours=rainbow(3)) + theme(axis.title.x = element_text(size=20)) + theme(axis.title.y = element_text(size=20)) + theme(axis.text.x = element_text(size=15)) + theme(axis.text.y = element_text(size=15)) + coord_cartesian(xlim=c(0.1,5000),ylim=c(100,2000000))

データに応じたグラフが得られるはず。このデータでは以下のようになった。最後に保存する。

ggsave(file="graph2.png")#セーブ

f:id:kazumaxneo:20170324134156j:plain