macでインフォマティクス

macでインフォマティクス

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

ゲノム情報はないが、モデル生物と近縁な生物のRNA seq 解析

 ゲノム情報はないが近縁種のゲノムが解読されているような生物でRNA seqを行うと決まったら、どんなワークフローで進めるべきだろうか?マイクロアレイと違い、RNA seqならde novoでも解析は不可能ではない。ゲノムがモデル生物とほぼ同じならば、深く考えなければ、なんだかできそうに思えてくる。

 実際、SRAのデータベースを見ているとそういったことを狙ったであろう解析データが無数に登録されている。その中にはpublish情報がないものもある。これが何を意味するかは置いておいて、モデル生物のゲノム情報の充実に伴い、そのゲノム情報を利用した近縁種の解析はこれからもますます増えてくると予想される。

ただしモデル生物に近縁でも塩基置換が起こってるだけでアライメントは難しくなる。ではde novo assemblyして進めるか?、というとせっかくゲノムが決まっているのだからなんとか利用したいという気持ちもあり、ジレンマに陥る。場合によりけりだと思うが、どういったワークフローが妥当なのか1つの例を通して考えてみることにする。

 

 今回はアラビドプシス属のA.kamchaticaのRNA seqデータを使う。シロイヌナズナとゲノムはかなり似ていると思われるが、ゲノムはまだ解読されていない。シロイヌナズナのecotypeでも塩基置換がかなり起きていることが知られており、A. kamchaticaのRNA seqデータがどれくらいナズナのリファレンスゲノムにアライメントされるかは未知数である。

 

今回はアセンブルしてリード数などの比較を行うが、それ以上は深入りせずプロトコルに抽象性を残しておく。その方が汎用性のあるプロトコルになると考える。

 

Arabidopsis. kamchatica (アラビドプシス属

f:id:kazumaxneo:20170730105553j:plain

 

 

データ

Transcriptome analysis of arabidopsis kamchatica under normal and stress conditions 

 

control1のデータ量

リード数 3600万 x 2

リード長 101-bp

ファイルサイズ 12.6GB x 2

 

 

万能なワークフローは存在しないが、今回はまずgenome guideのアセンブリと、de novoのアセンブリを実行し、それぞれの手法で作ったcDNA (contig)をTAIRのcDNAと比較する。比較するのは、生じたcDNAの数と長さ、core geneの存在数、TAIR cDNAとの相同性などになる。TAIR cDNAと比較してArabidopsis kamchatica固有の配列の割合を見積もることで、発現解析にはどのようなcDNAセットを使うのが妥当なのか分析する。(ここで言う優良なcDNA セットとは、単純に漏れが無くてダブりも無いセット)。

 

 

 

マッピングに使うアセンブルデータが決まったら、それをTAIR のcDNAとマージし、さらにcd-hit-estを使って95%以上同一の重複配列を廃除する。

 

 

 

 

 

1、contigの作成。

genomeガイドの方法1つと、de novo assemblyの手法2つをを心みる。

1、 genome-guided assembly

tophatでのmapping

bowtie2-build -f A_thaliana_TAIR10.fa chr
tophat2 -p 4 -I 50000 -i 20 -o tophat chr control1/DRR013381_1.fastq control1/DRR013381_2.fastq

 

マッピング率の分析

#アライメントされたリード
samtools flagstat accepted_hits.bam

15901037 + 0 in total (QC-passed reads + QC-failed reads)

1879156 + 0 secondary

0 + 0 supplementary

0 + 0 duplicates

15901037 + 0 mapped (100.00% : N/A)

14021881 + 0 paired in sequencing

7007096 + 0 read1

7014785 + 0 read2

6438422 + 0 properly paired (45.92% : N/A)

6950656 + 0 with itself and mate mapped

7071225 + 0 singletons (50.43% : N/A)

57762 + 0 with mate mapped to a different chr

22166 + 0 with mate mapped to a different chr (mapQ>=5)

 

#アライメントされなかったリード
samtools flagstat unmapped.bam

58401591 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 secondary

0 + 0 supplementary

0 + 0 duplicates

0 + 0 mapped (0.00% : N/A)

58401591 + 0 paired in sequencing

29204640 + 0 read1

29196951 + 0 read2

0 + 0 properly paired (0.00% : N/A)

0 + 0 with itself and mate mapped

0 + 0 singletons (0.00% : N/A)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

 ナズナゲノムにアライメントされたリード1590万、アライメントされなかったリード5840万。マッピング率はかなり悪い。この数値だとリファレンスベースで進めるとかなりの取りこぼしが出ると予想される。

 

ただし今回はアライメントの感受性は変えずこのまま進める。デフォルトのパラメータでどこまで妥当なデータが出せるか考える。

 

Cufflinks

cufflinks --no-update-check --overlap-radius 1 -o cufflinks accepted_hits.bam
  • –overlap-radius <int> Transfrags that are separated by less than this distance get merged together, and the gap is filled. Default: 50 (in bp).
  • –no-update-check Turns off the automatic routine that contacts the Cufflinks server to check for a more recent version.

 

Cuffmerge

echo cufflinks/transcripts.gtf > assemblies.txt
cuffmerge -s A_thaliana_TAIR10.fa assemblies.txt

/merged_asm にmerged.gtfができる。featureの数は61958だった。

ここからbedに変換してbedtoolsでfasta配列を抽出し、de novoで作ったfastaと比較する。

 

cat merged.gtf | awk '{OFS = "	"} {print $1,$4,$5,$3,$6,$7}' > merged.bed

#コマンド省略 perlスクリプトでbedの各featureのサイズを計算、7カラム目に追加= > sized.bed として出力

#サイズ順にソート
sort -k 7 -n -r sized.bed > size_sorted.bed

#コマンド省略 100bp以下を除去して保存

#bedtoolsでfastaを抽出。
bedtools getfasta -fi ../A_thaliana_TAIR10.fa -bed size_sorted.bed |fold -w 60 > cufflinks.fasta

featureの数は30422だった。大きく数が減ったのは、cuffmergeで0bpの配列が大量にできており、また100bp以下の非常に短いfeatureも大量にあったためである。

 

 

2、 de-novo assembly - rnaspades

rnaspades.py -t 20 -k auto -o RNA_spades -1 R1.fastq -2 R2.fastq

合計68374 contigできた。そのうち 100bp以上のcontigは68176だった。diploidのspeciesに使ってしまってますが、本来は推奨されていません、ご注意ください。

 

 

3、 de-novo assembly - trinity

Trinity --seqType fq --left R1.fastq --right R2.fastq --max_memory 40G --CPU 20

合計78657 contigできた。 そのうち 100bp以上のcontigは78656だった(1つだけ21 bp)。

アセンブルでエラーを起こした。フォーマットが古いせいだと思われたので、ヘッダーを修正して3行目を+に修正してランした。

 

4、 genome-guided assembly - trinity

trinityでも、ゲノムを参照しながらアセンブルする方法があるのを忘れていました。ここでは結果なしになりますが、リンクを載せて起きます。

https://github.com/trinityrnaseq/trinityrnaseq/wiki/Genome-Guided-Trinity-Transcriptome-Assembly

 

 

 

2、buscoによるcore geneの探索。

 

BUSCO -g input.fasta -o OUTPUT -l embryophyta_odb9 -m trans -c 8

1、 genome-guided assembly 

#BUSCO was run in mode: trans

 

Summarized benchmarks in BUSCO notation:

C:0%[D:0%],F:0%,M:0%,n:1440

 

13 Complete BUSCOs

13 Complete and single-copy BUSCOs

0 Complete and duplicated BUSCOs

113 Fragmented BUSCOs

1314 Missing BUSCOs

1440 Total BUSCO groups searched

13しかできていない。

 

2、 de-novo assembly - rnaspades

#BUSCO was run in mode: trans

 

Summarized benchmarks in BUSCO notation:

C:0%[D:0%],F:0%,M:0%,n:1440

 

937 Complete BUSCOs

668 Complete and single-copy BUSCOs

269 Complete and duplicated BUSCOs

164 Fragmented BUSCOs

339 Missing BUSCOs

1440 Total BUSCO groups searched

937見つかった。339 missingがある。

 

3、 de-novo assembly - trinity

#BUSCO was run in mode: trans

 

Summarized benchmarks in BUSCO notation:

C:0%[D:0%],F:0%,M:0%,n:1440

 

1078 Complete BUSCOs

797 Complete and single-copy BUSCOs

281 Complete and duplicated BUSCOs

105 Fragmented BUSCOs

257 Missing BUSCOs

1440 Total BUSCO groups searched

1078見つかった。257 missingがある。

 

参考; TAIR cDNA

Summarized benchmarks in BUSCO notation:

C:0%[D:0%],F:0%,M:0%,n:1440

 

1357 Complete BUSCOs

896 Complete and single-copy BUSCOs

461 Complete and duplicated BUSCOs

17 Fragmented BUSCOs

66 Missing BUSCOs

1440 Total BUSCO groups searched

1357見つかる。TAIRのcDNAにも66 missingがある。

 

まとめ

f:id:kazumaxneo:20170804143209j:plain

cuffmergeはかなり厳しい。このケースではゲノムガイドアセンブリは現実的ではないようだ。 trinityがベストのように思える。

 

 

#BUSCOは以前紹介しました。簡単な使い方をまとめています。

アセンブル結果をCore gene setの検出数で評価する BUSCO - macでインフォマティクス

 

 

3、contigの長さの比較

 

fastaのサイズについてヒストグラムを書き比較する。

check_fasta_length.pl RNA_spades/transcripts.fasta |cut -f 2 - > spades.txt #spadesのfastaの長さを抽出

check_fasta_length.pl trinity/transcripts.fasta |cut -f 2 - > trinity.txt #trinityのfastaの長さを抽出

check_fasta_length.pl Arabidopsis_thaliana.TAIR10.cdna.all.fa |cut -f 2 - > TAIRcDNA.txt #TAIRcDNAのfastaの長さを抽出

Rにて

spades=scan("spades.txt")
trinity=scan("trinity.txt")
cuffmerge=scan("cufflinks_contig")
TAIRcDNA=scan("TAIRcDNA.txt")

hist(spades, col = "#ff00ff40", breaks = 150)
hist(cuffmerge, col = "#00A47450", breaks = 150, add = TRUE) #グラフに追加
hist(TAIRcDNA, col = "#2F343720", breaks = 300, add = TRUE) #グラフに追加 

f:id:kazumaxneo:20170727011204j:plain

 黒がTAIRのcDNAの分布で、ピンクがrnaspadesでアセンブルしたcontigの分布、緑がcuffmergeでゲノムベースにアセンブルしたcontigの分布。縦軸は頻度。横軸は長さ(bp)。

cuffmergeで作られた領域は短いものが非常に多い。 

 

次にtrinityとrnaspadesの比較。

hist(spades, col = "#ff00ff40", breaks = 150) 
hist(trinity, col = "#0000ff40", breaks = 150, add = TRUE) #グラフに追加

 

f:id:kazumaxneo:20170804135811j:plain

ピンクがrnaspadesでアセンブルしたcontigの分布、青がTrinityでセンブルしたcontigの分布。Trinityの方が長いcontigが多い。

 

 

 

 

4、contigのblast解析

3つの方法で作成したcontigをTAIRのcDNAと比較する。blastnを使う。

 

まずblastのデータベースを作成。

makeblastdb -dbtype nucl -in Arabidopsis_thaliana.TAIR10.cdna.all.fa -out TAIR10_cDNA_db

 

相同性検索を実行。出力指定は可能だがdefaultで行う(リンク)。

#RNA_spades
blastn -query RNA_spades/transcripts.fasta -evalue 1e-10 -outfmt 6 -max_target_seqs 1 -db TAIR10_cDNA_db > RNA_spades_out.txt

#trinity
blastn -query trinity/transcripts.fasta -evalue 1e-10 -outfmt 6 -max_target_seqs 1 -db TAIR10_cDNA_db > trinity_out.txt
  • -outfmt 6 タブ区切り出力

出力されたファイルにはe-valueが1e-10以下を満たす配列が保存されている。中身を確認しておく。

1 2 3 4 5 6 7 8 9 10 11 12
AT1G19090.1 AT1G19090.1 99.88 2509 3 0 1 2509 1 2509 0 4617

1、2列目がクエリーとデータベースの名前である。さらに3列目が相同性(%)、4列目がhitした長さ、5列目がhitした領域中のミスマッチの数、6列目がギャップの数である。12列目にスコア (Bits) も出ている。

 

blastの結果をまとめる。48359あるTAIRのcDNAとヒットしたのは (e-value <1e-10)、

cuffmergeの30422 contigのうち30011 contig 

rnaspadesの68176 contigのうち50715 contig

Trinityの78657 contigのうち 63575 contig

 

hitしなかったcontigシロイヌナズナ ゲノムにはない配列だったか、シロイヌナズナゲノムに存在する配列だが、TAIR cDNAと転写領域が異なるためcDNAカタログ化されていない配列と考えられる。

 

 

 

 

 

2−4のまとめ

 ここまでの解析で、シロイヌナズナのゲノムベースでアセンブルするより、de novoでアセンブルするデータの方が漏れが少なく、cDNAのサイズも長くなることがわかった。

 ちなみに今回はリードの長さを分析してblast解析するという泥臭い方法をとったが、TransRateというツールを使うと、リファレンスとの相同性やエラー率など複数の側面からアセンブリを評価してくれる。

 

TransRateはこちらで紹介しています。

アセンブル結果をリードのアライメントパターンから評価する TransRate - macでインフォマティクス

  

 この後どうすべきかであるが、万能な方法はない。他の人はどうしているのだろう?例えばこの論文ではrudandantなcontigを除く方法を提案している。

 

 

 

最初書いた時はシロイヌナズナの情報を諦めきれず、contig配列とシロイヌナズナの配列を1つにしてcd-hit-estで選抜してより長い配列を使うことを考えていた。しかしよく考えてみると、cd-hitで選抜すると、完全長のcDNAになっているシロイヌナズナの配列が優先されてde novoで作ったcDNAは捨てられ、結局アライメントでうまく行かなくなってしまう可能性が高い。またそれこそ恣意的なデータセットを使っていると取られかねない。

 

 ということで、みんながやっているようにdenovo assmeblyで作ったcontigのみ使うことにする。解析に供した全サンプルのリードを使ってアセンブルすれば、完全なcDNAコレクションにはなつていないとしても、少なくともそのデータに関しては定量漏れの遺伝子が出る可能性は低い。

 

たたし、denovo assmeblyのデータだけ使うなら、アセンブルの精度が実験の品質を完全に左右することになるので、アセンブル精度を客観的に評価できる手法でアセンブル精度を担保しておく必要が出てくる。

 

以下は最近アクセプトされた論文での流れ。

 2016年 nature scientific report de nono transcriptomeのペーパー

1、trinityを使い、k-merが25、32でアセンブリ

2、できたcontigをBUSCOとTransRateで評価。

3、アセンブリには適切なk-mer coverageがあるとの先行研究から、ランダムサンプリングでリード数を減らし再アセンブル

4、1のデータと3のデータをBUSCOとTransRateで再評価

5、ランダムサンプリングしてリード数を減らしk-mer=32でアセンブルのがベストだと判断。

6、Trinotateでアノテーション

7、 blastで植物ゲノムとベストヒットしないものを排除。

 

 このように、k-mer、カバレッジ、などを試行錯誤してアセンブルし、それから評価ツールで良好なアセンブル条件を比較する流れを取っている。 

 

今回のケースで、近縁種のゲノムが出ていたとしても、直接鋳型としてgenome guide assemblyをするのは厳しいことを身をもって実感した。ゲノムガイドとde novoのどちらも検討してもよいが、最終的には、近縁種のモデル生物のゲノム情報はアノテーション情報として利用するのが精一杯と思われる。

 

 

追記1

下のリンク先の論文は、複数のアセンブルツールを併用し、cd-hit、TGI clustering tools、EvidentialGeneを使ってredundancyを下げることで、単一のアセンブルツールより高品質なRNAアセンブルを行なっている。

Combining Transcriptome Assemblies from Multiple De Novo Assemblers in the Allo-Tetraploid Plant Nicotiana benthamiana

こちらもそう

https://dl.sciencesocieties.org/publications/tpg/pdfs/8/1/plantgenome2014.10.0064

TGI clustering toolsのpaper

https://academic.oup.com/bioinformatics/article/19/5/651/239299

 

追記2

RNAのリファレンスガイドアセンブリで精度を上げる 

 

追記3

発現パラーンからクラスター化する。


追記4

アセンブルの評価

Transrate

BUSCO

 

 

追記5

Evaluation of de novo transcriptome assemblies from RNA-Seq data

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0553-5

 

追記6

Differential expression analysis of de novo assembled transcriptomes - Nadia Davidson from Australian Bioinformatics Network

 
追記7
ハイパフォーマンスなd novoアセンブラ


 

 

以下は、cd-hitを使ってシロイヌナズナのcDNAと統合しようと考えていた時の流れである。論文投稿のレベルを満たした手法とは思えないが、残しておく。 

 

5、マージ

ここから先は2つの手段を取る。

 

A、アセンブルして得られたcDNAからシロイヌナズナに見つからないcDNAだけ選り分け、それをTAIR cDNAとマージしてリファレンスとする。

B、アセンブルで得られたcDNAとTAIR cDNAをそのままマージし、CD-HIT-ESTで重複を除きリファレンスとする。

 

Aの方法を実行するには、シロイヌナズナのcDNAと相同性を示さないcDNAを集める必要がある。そこで、step4でblast解析した結果を相同性を示すcDNAのデータベースとして、no hitのcDNAをrnaspadesのアセンブルfastaから抽出する作業を行なった。そうゆうオプションやツールがありそうだが(絶対ある)、すぐには見つからなかったのでperlスクリプトを書いて処理した。

得られたno-hitのcDNA.fastaをblastにかける。organismはシロイヌナズナに限定し、最初の100 contigをblast解析する。e-value閾値はデフォルトでランした。

f:id:kazumaxneo:20170727144655j:plain

 

no-hitのcDNAだけをクエリにしているので、 スクリプトが上手くワークしているなら、シロイヌナズナのゲノムと全体にわたってマッチする配列は0のはずである。

 

結果は大まかに3パターンに分かれた。

 

highly homologous with A.thaliana genome

f:id:kazumaxneo:20170727150550j:plain

partially homologous with A.thaliana genome

f:id:kazumaxneo:20170727150611j:plain

no siginificant holology with A.thaliana genome

f:id:kazumaxneo:20170727150656j:plain

 予想通り多くのcDNAはナズナゲノムと高い相同性は示さなかった。このことは、Arabidopsis. kamchatica 固有の配列が抽出されていることを意味する。しかしながら、blastにかけた100 cDNAのうち1割ほどのcDNAはシロイヌナズナゲノムと配列全体に渡って高い相同性を示した。恐らくは、スプライシングバリアントというよりTAIRのcDNAカタログと全く異なる転写領域のため、blastで相同性を示さなかったのだと思われる。

シロイヌナズナのcDNAの再解析の論文も報告されているが、そういったデータと比べるとどうなるだろうか? http://www.biorxiv.org/content/biorxiv/early/2016/04/05/047308.full.pdf )

 

得られたno-hit CDAとTAIR 10 cDNAをマージして、no-hit_and_TAIRcDNA.faとした。step6のマッピングに進む。

 

 

次に Bの方法について書く。

まずrnaspadesで得られたfastaとTAIR cDNAを1つにして"TAIR_and_rnaspades_merged.fa"を作った。このファイルを入力としてcd-hit-estでクラスター化する。

cd-hit-est -i TAIR_and_rnaspades_merged.fa -o TAIR_and_rnaspades_merged_cd-hit_clustered.fa -c 0.95 -T 12
  • -i   input filename in fasta format, required
  • -o output filename, required
  • -c  sequence identity threshold, default 0.9 this is the default cd-hit's "global sequence identity" calculated as: number of identical amino acids in alignment
  • -T number of threads, default 1; with 0, all CPUs will be used.

閾値は95%とした。これ以上の相同性を示すcDNAは1つにまとめられる。

解析が終わると、クラスター化されて95%以上似た配列が除去されたfastaと、どのfastaクラスターになったかを示す~fa.clstrができる。中身を見てみる。

user$ head -10 cd-hit_est_clustered.fa.clstr 

>Cluster 2

0 15464nt, >AT3G02260.4... at +/98.96%

1 15700nt, >AT3G02260.2... *

2 15527nt, >AT3G02260.3... at +/99.98%

3 943nt, >NODE_18797_length_9... at -/98.09%

4 412nt, >NODE_35931_length_4... at +/99.03%

5 284nt, >NODE_45670_length_2... at +/97.18%

6 221nt, >NODE_54874_length_2... at -/96.38%

7 120nt, >NODE_65872_length_1... at -/99.17%

8 109nt, >NODE_67550_length_1... at +/93.58%

このクラスターだと8つのcDNAからなる。それがアスタリスクの ">AT3G02260.2 "にまとめられている。ただしcDNAの長さや塩基が編集を受けるわけではない。AT3G02260.2はクラスター化する前も後も15700-bpで100%同一の配列である。

 

cd-hit-estにより、116692あった配列が69650に減った。

ヒストグラムを書く。 

f:id:kazumaxneo:20170729112752j:plain

ピンクがcd-hit-estにかける前のcDNA、黒がcd-hit-estにかけた後のcDNAである。横軸はcDNAの長さ、縦軸が頻度である。特に短い配列が多く減っているように見える。

BUSCOを使い、 クラスター化でcore geneが減ってしまっていないか確認する。

BUSCO -g input.fasta -o OUTPUT -l embryophyta_odb9 -m trans -c 8

cd-hit-estにかける前のfasta

#BUSCO was run in mode: trans

 

Summarized benchmarks in BUSCO notation:

C:0%[D:0%],F:0%,M:0%,n:1440

 

1395 Complete BUSCOs

414 Complete and single-copy BUSCOs

981 Complete and duplicated BUSCOs

22 Fragmented BUSCOs

23 Missing BUSCOs

1440 Total BUSCO groups searched

TAIRのcDNAを含むので、core geneの保有数は1440のうち1395まで上がっている(rnaspades単独のデータだと937)。

 

cd-hit-estにかけた後のfasta

#BUSCO was run in mode: trans

 

Summarized benchmarks in BUSCO notation:

C:0%[D:0%],F:0%,M:0%,n:1440

 

1335 Complete BUSCOs

835 Complete and single-copy BUSCOs

500 Complete and duplicated BUSCOs

17 Fragmented BUSCOs

88 Missing BUSCOs

1440 Total BUSCO groups searched

core geneは60減って1335となった。ただしduplicateが大きく減り、その分、single copyが大きく増えた。これをどう捉えるか...

 

 

 

追記

この記事を書いた時点ではしりませんでしたが、非モデル生物のRNA seq解析なら、CORESTという1遺伝子由来のcontigをまとめてクラスター状態でカウント可能にするツールが使えそうです。


contigをまとめることの有用性は2016年のアライナーの論文でも指摘されています。"5 Application of quasi-mapping for clustering de novo assemblies" の段落を読んでみてください。

 

 

引用

SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing 

Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son Pham, Andrey D. Prjibelski, Alexey V. Pyshkin, Alexander V. Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev, and Pavel A. Pevzner.

Journal of Computational Biology. May 2012, 19(5): 455-477. 

 

TopHat: discovering splice junctions with RNA-Seq

Cole Trapnell Lior Pachter Steven L. Salzberg

Bioinformatics (2009) 25 (9): 1105-1111.

 

Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation

Cole Trapnell, Brian A Williams, Geo Pertea, Ali Mortazavi, Gordon Kwan, Marijke J van Baren, Steven L Salzberg, Barbara J Wold & Lior Pachter Nature

Biotechnology 28, 511–515 (2010) doi:10.1038/nbt.1621

 

 

The Sequence Alignment/Map format and SAMtools.

Li H1, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup.

Bioinformatics. 2009 Aug 15;25(16):2078-9. doi: 10.1093/bioinformatics/btp352. Epub 2009 Jun 8.

 

CD-HIT Suite: a web server for clustering and comparing biological sequences.

Ying Huang, Beifang Niu, Ying Gao, Limin Fu and Weizhong Li.

Bioinformatics, 2010(26): 680-682.full text

 

Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences.

Weizhong Li and Adam Godzik.

Bioinformatics, 2006(22): 1658-1659. full text

 

Tolerating some redundancy significantly speeds up clustering of large protein databases.

Weizhong Li, Lukasz Jaroszewski and Adam Godzik.

Bioinformatics, 2002(18): 77-82. full text

 

Clustering of highly homologous sequences to reduce the size of large protein databases.

Weizhong Li, Lukasz Jaroszewski and Adam Godzik.  

Bioinformatics, 2001(17): 282-283. full text

 

 

 

フォーマット変換 FASTA => BED

 

fasta.faiから作る。

samtools faidx input.fasta 
awk '{print $1 "\t0\t" $2}' input.fasta.fai > output.bed

 

またはpythonスクリプトを使う。

pip install pyfaidx 
faidx --transform bed input.fasta > output.bed

 

ヒトゲノムhg19ならこのようなbedができる。

user$ head hg19.bed 

chr1 0 249250621

chr2 0 243199373

chr3 0 198022430

chr4 0 191154276

chr5 0 180915260

chr6 0 171115067

chr7 0 159138663

chr8 0 146364022

chr9 0 141213431

chr10 0 135534747

 

 

 

引用

Shirley MD, Ma Z, Pedersen B, Wheelan S. Efficient "pythonic" access to FASTA files using pyfaidx. PeerJ PrePrints 3:e1196. 2015.

https://github.com/mdshw5/pyfaidx

 

Biostarスレッド

https://www.biostars.org/p/191052/

BEDフォーマット

 

UCSCのゲノムブラウザーなどで使うフォーマット。最初の3列が必須で、オプションでさらに9列情報がつく場合がある。

 

最初の3列に記載する情報

  1. クロモソームの名前(e.g., chr1)
  2. リードや遺伝子のスタートポジション(ポジションは1でなく0スタート
  3. リードや遺伝子のエンドポジション

 

追加の9列

  1. 名前
  2. 1-1000のスコア。スコアに応じて以下のようなカラー情報を持たせる事も可能である。f:id:kazumaxneo:20170703222731j:plain
  3. リードや遺伝子の向き(+/-)
  4. CDSのスタートポジション。リードなら2の座標と同じになる。
  5. CDSのエンドポジション。
  6. exonの数
  7. 各exonのサイズ(数値をコンマで区切り全て記載する)
  8. exonのスタート位置。

  

 BAMから変換したり、bed同士のオーバラップなどを分析するにはbedtoolsを使います。bedtoolsは以下で紹介しています。


引用

フォーマット一覧(UCSC

https://genome.ucsc.edu/FAQ/FAQformat.html

アメリエフのブログ(旧)

http://blog.amelieff.jp/?eid=195350

 

 

 

VCF, GTF, GFF などを BED に変換 する BEDOPS

2019 6/17 追記

2020 2/21 タイトル修正

2020 3/30 help追記

 

BEDヘの変換はawkperlpythonスクリプトで簡単にできるが、BEDOPSのvcf2nedを使うと、indelの種類などによってフィルタリングしながら分類することができ便利である。

 

インストール

#homebrew
brew install BEDOPS

#bioconda(link)
conda install -c bioconda -y bedops

bedops

$ bedops 

bedops

  citation: http://bioinformatics.oxfordjournals.org/content/28/14/1919.abstract

            https://doi.org/10.1093/bioinformatics/bts277

  version:  2.4.37 (typical)

  authors:  Shane Neph & Scott Kuehn

 

      USAGE: bedops [process-flags] <operation> <File(s)>*

 

          Every input file must be sorted per the sort-bed utility.

          Each operation requires a minimum number of files as shown below.

            There is no fixed maximum number of files that may be used.

          Input files must have at least the first 3 columns of the BED specification.

          The program accepts BED and Starch file formats.

          May use '-' for a file to indicate reading from standard input (BED format only).

 

      Process Flags:

          --chrom <chromosome> Jump to and process data for given <chromosome> only.

          --ec                 Error check input files (slower).

          --header             Accept headers (VCF, GFF, SAM, BED, WIG) in any input file.

          --help               Print this message and exit successfully.

          --help-<operation>   Detailed help on <operation>.

                                 An example is --help-c or --help-complement

          --range L:R          Add 'L' bp to all start coordinates and 'R' bp to end

                                 coordinates. Either value may be + or - to grow or

                                 shrink regions.  With the -e/-n operations, the first

                                 (reference) file is not padded, unlike all other files.

          --range S            Pad or shrink input file(s) coordinates symmetrically by S.

                                 This is shorthand for: --range -S:S.

          --version            Print program information.

 

      Operations: (choose one of)

          -c, --complement [-L] File1 [File]*

          -d, --difference ReferenceFile File2 [File]*

          -e, --element-of [bp | percentage] ReferenceFile File2 [File]*

                 by default, -e 100% is used.  'bedops -e 1' is also popular.

          -i, --intersect File1 File2 [File]*

          -m, --merge File1 [File]*

          -n, --not-element-of [bp | percentage] ReferenceFile File2 [File]*

                 by default, -n 100% is used.  'bedops -n 1' is also popular.

          -p, --partition File1 [File]*

          -s, --symmdiff File1 File2 [File]*

          -u, --everything File1 [File]*

          -w, --chop [bp] [--stagger <nt>] [-x] File1 [File]*

                 by default, -w 1 is used with no staggering.

      

Example: bedops --range 10 -u file1.bed

      NOTE: Only operations -e|n|u preserve all columns (no flattening)

 

 

 

公式マニュアル

http://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html

 

ラン

vcfからbedに変換する。

vcf2bed < gatk.vcf > gatk.bed
  • --do-not-sort (-d)  Do not sort BED output with sort-bed
  • --snvs (-v) Report only single nucleotide variants
  • --insertions (-t) Report only insertion variants
  • --deletions (-n) Report only deletion variants
  • --keep-header (-k) Preserve header section as pseudo-BED elements

snpsのみbedに変換する。

vcf2bed --snvs < gatk.vcf > gatk_snps.bed

 

塩基置換、挿入、欠損の数を数える。

vcf2bed --snvs < gatk.vcf|wc -l       #SNV
vcf2bed --insertions < gatk.vcf|wc -l  #Insertion
vcf2bed --deletions < gatk.vcf|wc -l  #Deletion

  

vcf2bedはBAM、GFF、GTF、GVF、PSL、RepeatMasker (OUT)、SAM、VCF、WIGなど多様なフォーマットをBEDに変換することができる。

GFF(GFF3)をbedに変換する。

convert2bed --input=gff < input.gff3 > output.bed

 

 またはawkを使う。以下のようにして6列フォーマットのBEDに変換できる。

cat input.gtf | awk '{OFS = "\t"} {print $1,$4,$5,$3,$6,$7}' > output.bed

 awkはデフォルトスペース区切り出力だが、bedtoolsはタブを区切りとして認識するので、タブ区切りを指定。

 

 

追記

BEDからGTF (cent OSで動作確認)

awk '{print $1"\t"$7"\t"$8"\t"($2+1)"\t"$3"\t"$5"\t"$6"\t"$9"\t"(substr($0, index($0,$10)))}' input.bed > output.gtf

 

 

 

BEDを使って何かするにはbedtoolsを使います。


 

 

引用

BEDOPS: high-performance genomic feature operations
Neph S1, Kuehn MS, Reynolds AP, Haugen E, Thurman RE, Johnson AK, Rynes E, Maurano MT, Vierstra J, Thomas S, Sandstrom R, Humbert R, Stamatoyannopoulos JA.

Bioinformatics. 2012 Jul 15;28(14):1919-20 

 

How To Convert Bed Format To Gtf?

 

bedtoolsでpromoter配列を抽出する。

 

Biostarとbedtoolsの公式サイトに、bedtoolsを使ったpromoter配列の抽出の仕方がまとめられている(How To Use Bedtools To Extract Promoters From A Mouse Bed File)。試してみる。

 

このようなbedファイルの各featureの上流を抽出する。

user$ head -5 genes.bed 

chr1 11873 12227 NR_046018_exon_0_0_chr1_11874_f

chr1 12612 12721 NR_046018_exon_1_0_chr1_12613_f

chr1 13220 14409 NR_046018_exon_2_0_chr1_13221_f

chr1 14361 14829 NR_024540_exon_0_0_chr1_14362_r

chr1 14969 15038 NR_024540_exon_1_0_chr1_14970_r

 

 

まずクロモソームのサイズを書いたファイルを準備する。単純なタブ区切りのテキストファイルであればよい。ここではsamtools faidxを使いfastaから作っている。

samtools faidx reference.fa
cut -f 1,2 reference.fa.fai > genome.txt

 

上流領域を抽出する。bedtoolsで隣接 featureを取れるコマンドは複数あるが、今回はflankコマンドを使う。

bedtools flank -i genes.bed -g genome.txt -l 1000 -r 0 -s > genes.1kb.promoters.bed

flankコマンドで隣接領域を抽出している。この時-l 1000 -r 0とすることで隣接する上流1000-bp、下流0-bpが抽出されるように工夫している。ちなみにクロモソームの末端でクロモソームからはみ出すfeatureがある場合、トリミングして出力される。

 

さらにbedtoolsのgetfastaで配列を抽出する。

bedtools getfasta -fi reference.fa -bed genes.1kb.promoters.bed |fold -w 60 > genes.1kb.promoters.bed.fa
  • -fi Input FASTA file
  • -bed BED/GFF/VCF file of ranges to extract from -fi

-fiでリファレンスゲノムを指定して使う。また、パイプを使いfoldで60文字ごとに改行している。少し重くなるので、不要なら外す。上記コマンドはヘッダー名が違うとワークしないので気をつける。

出力は、例えば grep -n "^>" input.fasta | wc -l で配列数をカウントできる。

 

 

使用したpromoter配列の塩基組成や偏りなどの基本情報をまとめたい場合、bedtoolsのnucを使う。

先ほど作ったbedを使い、以下のように打つ。

bedtools nuc -fi reference.fa -bed genes.1kb.promoters.bed > ATGC.bed

カラムの後半に以下のような情報がつく。

  • 1) %AT content
  • 2) %GC content
  • 3) Number of As observed
  • 4) Number of Cs observed
  • 5) Number of Gs observed
  • 6) Number of Ts observed
  • 7) Number of Ns observed
  • 8) Number of other bases observed
  • 9) The length of the explored sequence/interval.
  • 10) The seq. extracted from the FASTA file. (opt., if -seq is used)
  • 11) The number of times a user's pattern was observed. 

  

平均や最大/最小などの情報はRでまとめる。

 

 

上で使ったbedtoolsのコマンドはこちらでまとめています。

 

 

引用

How To Use Bedtools To Extract Promoters From A Mouse Bed File

オオムギ(大麦)のRNA seq解析

 勉強会用資料

時系列データ

 

publishされた論文

http://www.sciencedirect.com/science/article/pii/S1631069115000888?via=ihub

 

利用するシーケンスデータ

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP032854

 

fastaとgtf

http://plants.ensembl.org/Hordeum_vulgare/Info/Index

 

シーケンスデータのダウンロード(上記リンクよりAccession Listをダウンロードして使用)

prefetch --option-file SRR_Acc_List.txt  --max-size 1000000000

#SRAのID、例えばSRR6364639を直接ダウンロードするなら
prefetch SRR6364639 --max-size 1000000000
  •   --max-size 1000000000 20GB以上のファイルのダウンロードに必要。

 

fastqに変換 

fastq-dump <inputt.sra> -O <directory>
  • -O 指定したディレクトリに出力。
  • --split-files  paired-endならをつける。

 

 index

bowtie2-build -f Hordeum_vulgare.fa chr

 

mapping

tophat2 -I 500000 -i 20 -p 15 --read-mismatches 5 --read-gap-length 4 --max-multihits 100 --read-edit-dist 6 -G Hordeum_vulgare.Hv_IBSC_PGSB_v2.36.gtf -o tophat_Bs1_2 chr Bs3_24/SRR1049655.fastq

 

featurecount

featureCounts -T 8 -t exon -g gene_id -a Hordeum_vulgare.Hv_IBSC_PGSB_v2.36.gtf -o geaturecounts_all.txt tophat_Bs1_2/accepted_hits.bam tophat_Bs2_2/accepted_hits.bam tophat_Bs3_2/accepted_hits.bam tophat_Bs1_24/accepted_hits.bam tophat_Bs2_24/accepted_hits.bam tophat_Bs3_24/accepted_hits.bam 

 

 

 raw read

library("ggplot2") 
library("reshape2")
rawCount = read.table("featureCounts_counts_trimmed.txt", header = TRUE, row.names = 1)
head(rawCount) #必ず毎回確認
artificialCount = log2(rawCount + 1)
head(artificialCount)
df <- melt(logdata) #data.frame形式に変換。
head(df) #必ず毎回確認
g <- ggplot(df, aes (x = variable, y = value )) + geom_boxplot()
plot(g)

 

f:id:kazumaxneo:20170716161733j:plain

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

bedtools

追記 bedgraph出力 

2019 9/4 インストール、twitterリンク追加

 

BEDファイルのオーバーラップ領域を抽出したり、マージしたりできるツール。BED以外にGFF、VCFも扱うことができる。bedtools <command> -a .bed -b .bedという使い方が基本。-aで指定したbedを-bで指定したbedと比較する。出力はリダイレクト(>)で指定しないと標準出力となる。中にはGC/AT含量を調べるコマンドまであり、多岐に渡っている。使えそうなコマンドを紹介してみます。

f:id:kazumaxneo:20180324100345j:plain

公式より転載。

 

 

 

インストール 

#bioconda (link)
mamba install -c bioconda -y bedtools

 

コマンド

intersect

オーバーラップする領域の抽出。

aのfeatureのうちbの遺伝子とオーバーラップする領域だけをコールする。

bedtools intersect -a reads.bed -b genes.bed -sorted > output.bed
  • -wa Write the original entry in A for each overlap.
  • -wo Write the original A and B entries plus the number of base pairs of overlap between the two features.
  • -c For each entry in A, report the number of overlaps with B.
  • -f Minimum overlap required as a fraction of A. - Default is 1E-9 (i.e., 1bp). - FLOAT (e.g. 0.50)

-waをつけるとbとオーバーラップするaの領域がそのまま出力され、-waをつけないとaとbでオーバーラップする領域が出力される(参考)。-woをつけると右端カラムにオーバーラップする領域のサイズ付きで出力される。-cをつけると右端カラムにオーバーラップする領域数付きで出力される。デフォルトでは1bp以上のオーバーラップを重複とみなすが、-fで指定することでその閾値をコントロールできる。厳しくするなら、例えば-f 0.5にすれば50%以上オーバーラップを重複に設定できる。

 

(複数ファイル同時)オーバーラップする領域の抽出。

aのfeatureのうち3つのファイルのオーバーラップする領域をコールする。

bedtools intersect -a reads.bed -b genes1.bed genes2.bed genes3.bed -sorted > output.bed

-wa -wb をつけると重複したfeature数がわかる。-namesをつけると、どのデータ由来かわかる(-names gene1 gene2 gene3 )。

 

オーバーラップしない領域の抽出。

aのfeatureのうちbの遺伝子領域にない領域だけをコールする。

bedtools intersect -a reads.bed -b genes.bed -v -sorted > output.bed

 

Aとオーバーラップし、Bとオーバーラップしない領域の抽出。

aのfeatureのうちLINEとオーバーラップしており、かつSINEとはオーバーラップしていないfeatureを抽出。

bedtools intersect -a genes.bed -b LINES.bed |  
bedtools intersect -a stdin -b SINEs.bed -v > output.bed

パイプと標準入力(stdin)を使い2回連続比較している。順次実行してもよい。

  

merge

オーバーラップしているfeatureを統合して返す。

bedtools merge -i reads.bed > output.bed
  • -d Maximum distance between features allowed for features to be merged. - Def. 0. That is, overlapping & book-ended features are merged.

-dはdefaultでは0だが、1bp以上に指定することでオーバーラップしてないが-d指定値以内にあるfeatureを統合できる(-d 1000なら1000bp以内のfeatureを統合できる)。 -c 4 -o count,collapseをつければ、どのfeatureがどのfeatureと何回オーバーラップしたのかも出力される。

 

マージコマンドは、chr順、スタートポジション順にソートされていないとエラーを起こす。unixのsortコマンドでソートしてから使う(sort -k1,1 -k2,2n リンク)。似たコマンドにclusterがある。clusterはマージを実行せず、そのfeatureについてタグIDをつけて返す。

 

complement

featureされている領域以外の部位を返す。exon配列を入力すれば、intronとintergenicの領域が出力される。使うにはgenome配列情報が必要。

bedtools complement -i exon.bed -g genome.txt > non-exonic.bed

genomeはタブで仕切られた以下のような長さが記載してあるファイル。

User$ head -5 genome.txt 

chr1 249250621

chr10 135534747

chr11 135006516

chr11_gl000202_random 40103

chr12 133851895

samtools faidxで簡単に作れる。

samtools faidx reference.fa
cut -f 1,2 reference.fa.fai > genome.txt

complementと似たコマンドにsubstractがある。

 

substract

featureされている領域以外の部位を返す。exon配列を入力すれば、intronとintergenicの領域が出力される。complementコマンドと異なり、bはbedファイルを指定する。

bedtools complement -a genes.bed -b introns.bed > exons.bed

例えばintron (/exon) を-bで指定して、exon (/intron) を抽出するような用途に使う。

 

flank

featureされている領域に隣接した領域を返す。

bedtools flank -i gens.bed -g genome.txt -b 10 > intergenic.bed

-bで抽出する配列の長さを指定する。-b 10ならfeatureの上流、下流10-bpずつが出力される。-l 2 -r 3ならfeatureの上流2-bp、下流3-bpが出力される。

 

slop

featureされている領域を指定値だけ増減させて返す。

bedtools slop -i reads.bed -g genome.txt -b 5 > output.bed
  • -b Increase the BED/GFF/VCF entry -b base pairs in each direction. - (Integer) or (Float, e.g. 0.1) if used with -pct.
  • -l The number of base pairs to subtract from the start coordinate. - (Integer) or (Float, e.g. 0.1) if used with -pct.
  • -r The number of base pairs to add to the end coordinate. - (Integer) or (Float, e.g. 0.1) if used with -pct.
  • -s Define -l and -r based on strand. E.g. if used, -l 500 for a negative-stranded feature, it will add 500 bp downstream. Default = false.
  • -pct Define -l and -r as a fraction of the feature's length. E.g. if used on a 1000bp feature, -l 0.50, will add 500 bp "upstream". Default = false.

-bで抽出する配列の長さを指定する。-b 5ならfeatureの上流末端と下流末端がそれぞれ5-bpずつ伸びて出力される。”5 100”なら”0 105”になる。-l 2 -r 3ならfeatureの上流2-bp、下流3-bpずつ伸びて出力される。-sもつけると、-l と-rで伸びる領域がリードの方向性を考慮して決められる。

 

genomecov

genome全体に対するfeatureのカバレッジを返す。mergeコマンドと同様、bedはあらかじめソートされている必要がある。"-i"でbedを指定する場合はgenome配列情報も必要。"-ibam"でbamを指定する場合、genome配列情報は不要。

bedtools genomecov -i exons.bed -g genome.txt > intergenic.bed

#0カバレッジの領域だけ出力(awkフィールド4が1以下か判定)
bedtools genomecov -bga -ibam sort.bam | awk '$4 < 1' > zero-coverage
#0カバレッジの領域だけ出力。awkの代わりにgrepを使う
bedtools genomecov -bga -ibam sort.bam | grep -w 0$ > zero-coverage

このコマンドは他より時間がかかる。-bgをつけるとfeatureのカバレッジを返す。-bgをつけないとゲノム全体のカバレッジを一定の間隔で区切って計算した値が返される。

  • -ibam  The input file is in BAM format. Note: BAM _must_ be sorted by position
  • -bg Report depth in BedGraph format. For details, see: genome.ucsc.edu/goldenPath/help/bedgraph.html
  • -bga   Report depth in BedGraph format, as above (-bg). However with this option, regions with zero coverage are also reported. This allows one to quickly extract all regions of a genome with 0 coverage by applying: "grep -w 0$" to the output.
  • -strand   Calculate coverage of intervals from a specific strand. With BED files, requires at least 6 columns (strand is column 6). - (STRING): can be + or -

  • -pc   Calculate coverage of pair-end fragments. Works for BAM files only

 似た名前のコマンドにcoverageがあるが、coverageコマンドはbed同士を比較して被っている割合、個数などを出す。

bedgraph出力

bedtools genomecov -i exons.bed -ibam input.bam -bg > bedgraph_output

 

coverage

aのfeatureについてbのfeatureとオーバーラップしている領域の割合、個数などの情報を返す。

bedtools coverage -a A.bed -b B.bed > coverage.bed

公式マニュアルの図がわかりやすいのでそのまま引用する。

Chromosome  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

BED FILE A  ***************     ***************     ******    **************

BED File B  ^^^^ ^^^^              ^^             ^^^^^^^^^    ^^^ ^^ ^^^^
              ^^^^^^^^                                      ^^^^^ ^^^^^ ^^

Result      [  N=3, 10/15 ]     [  N=1, 2/15 ]     [N=1,6/6]   [N=6, 12/14 ]

上のaのfeatureについて、bのfeatureを使いcoverageを分析した結果はresultのようになる。-sをつけるとfeature同士の向きが同じものだけカウントされる。

 

closest

aのそれぞれのfeatureについて、最も近いbのfeatureをコール。

bedtools closest -a genes.bed -b ALUs.bed > output.bed

  -t last-t firstなどをつけると、同じ距離に複数featureがあった場合の選択をコントロールできる。

  • -t How ties for closest feature are handled. This occurs when two features in B have exactly the same "closeness" with A. By default, all such features in B are reported. Here are all the options:

   - "all"    Report all ties (default).

   - "first" Report the first tie that occurred in the B file. -

   - "last" Report the last tie that occurred in the B file.

 

 

bamtobed

 bamをbedに変換する(bed)。 

bedtools bamtobed -i reads.bam > reads.bed

bamをbedpeに変換する(bedpe)。 

bedtools bamtobed -i reads.bam -bedpe > reads.bedpe

似たコマンドにBEDOPSのgtf2bedがある。gtf2bedの方がより様々なフォーマットからBEDを作成可能である。

フォーマット変換 GTF => BED - macでインフォマティクス

 

bamtofastq

bamをfastqに変換する。

bedtools bamtofastq -i reads.bam -fq single.fastq
  • -fq2 FASTQ for second end. Used if BAM contains paired-end data. BAM should be sorted by query name (samtools sort -n aln.bam aln.qsort) if creating paired FASTQ with this option.

 ペアリードのbamなら-fq2で相方の出力も指定する。

 

bedtobam

bedをbamに変換する。

bedToBam [OPTIONS] -i reads.bed -g <GENOME> > reads.bam
  • -mapq Set a mapping quality (SAM MAPQ field) value for all BED entries. Default: 255

似たコマンドにbedpetobamリンク)がある。

 

getfasta

bedからfasta配列を出力する。 

bedtools getfasta -fi ref.fasta -bed reads.bed > reads.fasta
  • -fi Input FASTA file
  • -s Force strandedness. If the feature occupies the antisense strand, the sequence will be reverse complemented. Default: strand information is ignored.
  • -bedOut Report extract sequences in a tab-delimited BED format instead of in FASTA format.
  • -name Use the “name” column in the BED file for the FASTA headers in the output FASTA file.
  • -fo output file (リダイレクトでも同じ)

bedで指定した領域の配列が出力される。配列はfeatureの向きを考慮して出力されるので、リファンレスと同じ配列の強制出力には-sをつける。

コマンドを使う場合、リファレンスのfastaのヘッダーはbedと同じものである必要がある。また、出力されるfastaは改行されない。長くて見にくい時はfoldコマンドで折り返す。( | fold -w 60 > reads.fasta)

 

window

CNVの上流/下流10000bp以内にあるgeneを調べる。

bedtools window -a CNVs.bed -b genes.bed -w 10000

-l -rをつけることで、上流と下流の数値を個別に指定できる。例えば-l 10000 -r 5000なら上流10000bp以内、下流5000bp以内となる。defaultは1000bp。smとSmをつけると向きも考慮できる。

  • -sm Only report hits in B that overlap A on the _same_ strand. - By default, overlaps are reported without respect to strand.
  • -Sm Only report hits in B that overlap A on the _opposite_ strand. - By default, overlaps are reported without respect to strand.

 

annotate

オーバーラップしているfeatureの個数、名前などの情報を返す。

bedtools annotate -both -i variants.bed -files genes.bed conserve.bed known_var.bed
  • -both Report the counts followed by the % coverage. - Default is to report the fraction of -i covered by each file.
  • -s Require same strandedness. That is, only counts overlaps on the _same_ strand. - By default, overlaps are counted without respect to strand.
  • -S Require different strandedness. That is, only count overlaps on the _opposite_ strand. - By default, overlaps are counted without respect to strand.

例えば変異をコールしたファイルのfeatureをgene.bedと比較して、コールされた部位のアノテーションをつけることに使える。上の例では3つのbedと比較し、同時に3つのアノテーションをつけている。-sや-Sを指定することでリードの向きを限定できる。

 

groupby 

featureをグループ分けして、様々な情報を返す。

bedtools groupby [OPTIONS] -i <input> -g <group columns> -c <op. column> -o <operation>

タブ仕分けされたファイルならなんでも使えるため、 公式ページではbed以外のテキストファイルの分析ツールとしての利用価値もあるとしている。

 作成途中

 

maskfasta

指定した領域以外をNにしてfastaを返す。

bedtools maskfasta -fi ref.fasta -bed genes.bed -fo output.fasta
  •  -fi  Input FASTA file
  • -bed BED/GFF/VCF file of ranges to mask in-fi
  • -fo    Output FASTA file
  • -soft Enforce "soft" masking. That is, instead of masking with Ns, mask with lower-case bases.

 -softをつけると、指定した領域以外をNでなく小文字にして返す。

 

shift

featureを指定した数値だけ移動させる。

bedtools shift -i reads.bed -g genome.txt -s 5 > reads_shifted.bed
  • -s Shift the BED/GFF/VCF entry -s base pairs. Integer or Float (e.g. 0.1) if used with -pct.
  • -m Shift entries on the - strand -m base pairs. Integer or Float (e.g. 0.1) if used with -pct
  • -p Shift entries on the + strand -p base pairs. Integer or Float (e.g. 0.1) if used with -pct
  • -pct Define -s, -m and -p as a fraction of the feature’s length. E.g. if used on a 1000bp feature, -s 0.50, will shift the feature 500 bp “upstream”. Default = false.

-s 5なら5bpポジションが下流にシフトする。-p 2 -m -3なら+のfeatureは2-bp下流、-のfeatureは3-bp上流にそれぞれシフトする。-s 0.5 -pctなら100-bpのfeatureが50-bp下流にシフトする。

  

nuc

featureの塩基組成を分析して返す。

bedtools nuc -fi input.fasta -bed genes.bed > oputput.bed

以下の情報が出力される。

  • 1) %AT content
  • 2) %GC content
  • 3) Number of As observed
  • 4) Number of Cs observed
  • 5) Number of Gs observed
  • 6) Number of Ts observed
  • 7) Number of Ns observed
  • 8) Number of other bases observed
  • 9) The length of the explored sequence/interval.
  • 10) The seq. extracted from the FASTA file. (opt., if -seq is used)
  • 11) The number of times a user's pattern was observed. 

GC/AT含量の分析に広く活用できると思われる。

 

random

指定された配列からrandomにfeatureを発生させる。

bedtools random -g genome.txt |sort -k1,1 -k2,2n > randam.bed
  • -l The length of the intervals to generate. Default = 100
  • -n The number of intervals to generate. Default = 1,000,000

部位がランダムに選ばれるので、sortコマンドに渡し、名前とポジションでソートして出力した方が見やすくなる(sort -k1,1 -k2,2n リンク。 

デフォルトではfeature間の間隔が 1,000,000なので、ゲノムが小さいとたくさんfeatureが出来すぎてしまう。-nを使いゲノムサイズに合わせて調節する(featureのサイズは-lで変えられるが、サイズをランダムにすることは出来ないらしい。-nで少しオーバーラップするくらいにfeatureを発生させて、bedtools merge -dすれば少しランダムなサイズのfeatureになる。しかし真のランダムとは程遠い...)。

 

 

 他にも様々なコマンドがあります。fisher検定して有意なのか偶然なのか判定させたり、オーバーラップの度合いを計算したりするなど多種多様なコマンドがあります。興味があれば公式ページから探してみてください。

The BEDTools suite — bedtools 2.26.0 documentation

 

その他、公式ページではbedtoolsを使った応用例がいくつか記載されています。exome解析やRNA seq解析のカバレッジを調べる方法などを知ってると便利です。

http://bedtools.readthedocs.io/en/latest/

時間があればここでもまとめてみます。

 

 

 

BEDフォーマットは別に紹介しています。

 

BEDがない場合、GTFやGFFから変換して作ります。

 

追記

bedutils


こちらも確認してください(主にヒト)。


 

 

引用

BEDTools: a flexible suite of utilities for comparing genomic features.

Quinlan AR1, Hall IM.

Bioinformatics. 2010 Mar 15;26(6):841-2. 

https://www.ncbi.nlm.nih.gov/pubmed/20110278

 

Example usage

http://bedtools.readthedocs.io/en/latest/content/example-usage.html