題:
如何將FASTA轉換為BED
SmallChess
2017-05-18 05:14:58 UTC
view on stackexchange narkive permalink

我有一個FASTA文件：

  > Sequence_1GCAATGCAAGGAAGAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC ... > Sequence_2GCAATGCAAGGAAGAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCpre.A。每個序列的BED文件如下： 
  Sequence_1 0 1500Sequence_2 0 1700  
 
 BED區域將僅僅是序列的大小。
 
 問：我之前是通過單行命令執行的。我不記得那是什麼，那是在Biostars上。我現在找不到帖子。進行轉換的最簡單方法是什麼？
你是說這個生物之星的帖子嗎？ https://www.biostars.org/p/15476/
@Greg帖子沒有給我一行命令來執行我想要的操作。我100％肯定我前一段時間在生物明星的某個地方看到過它，但是我不記得可以做到這一點的程序的名稱。
這個如何？ https://www.biostars.org/p/15476/#189089
@Greg這是一個Python程序。我不記得這個名字，因此也不知道如何搜索。
我們是否需要創建一個新標籤，例如“格式轉換”，“文件轉換”？如我所料，將會有更多的問題從“ X轉換為Y格式”類型。
七 答案:
bli
2017-05-18 16:33:24 UTC
view on stackexchange narkive permalink

您可以使用 bioawk輕鬆地做到這一點，它是awk的一個版本，具有促進生物信息學的附加功能：

  bioawk -c fastx'{print $ name“ \ t0 \ t“ length（$ seq）}'test.fa

-c fastx 告訴程序應該將數據解析為fasta或fastq格式。這使得awk命令中的 $ name 和 $ seq 變量可用。

Sam Nicholls
2017-05-18 13:40:08 UTC
view on stackexchange narkive permalink

將 FASTA 編入索引是一種很好的做法，因此您可以利用可能已經擁有的 .fai 。如果不是，您可以使用 samtools 生成索引，並使用一些 awk 製作您的 BED ：

  samtools faidx $ fastaawk'開始{FS =“ \ t”}; {print $ 1 FS“ 0” FS $ 2}'$ fasta.fai > $ fasta.bed

這將保持製表符分隔，但是您可以刪除 BEGIN 語句使用空格。 BED規範僅需要簡單的 BED 格式的“空格”。

neilfws
2017-05-18 06:01:14 UTC
view on stackexchange narkive permalink

您可以修改此單線襯套。請注意，它假定序列ID的長度不超過100個字符，並且在標題行的序列ID後面沒有描述。

  cat myseqs.fasta | awk'$ 0〜“ >” {打印c; c = 0; printf substr（$ 0,2,100）“ \ t0 \ t”; } $ 0！〜“ >” {c + = length（$ 0）;} END {打印c; }'

否則，任何Bio *庫（Perl，Python，Ruby）都提供FASTA格式的解析器，該解析器將提取序列ID和長度。

我要指出嚴格說來，這與BED類似，但並不是完全如此，因為BED映射到染色體或更長序列對像上的坐標。

Scott Gigante
2017-05-18 06:08:32 UTC
view on stackexchange narkive permalink

受此答案的啟發，對一個有關讀取長度分佈的相關問題的啟發，您可以使用Biopython做到這一點：

 來自Bio.SeqIO import parsewith open（“ regions .bed“，” w“）作為bed：用於parse（” regions.fasta“，” fasta“）中的記錄：print（record.id，0，len（record.seq），sep =” \ t“，文件= bed）
我喜歡Python。我更喜歡這個答案，因為它是在Python中使用的，這就是我們在2017年應該使用的答案。但是為什麼這看起來像Perl？可讀性很重要。我想提出一個更簡潔的代碼段，但是看起來註釋不允許代碼。編輯您的答案會被認為是一種很好的語氣嗎？
點了！我已經試過了，但是歡迎您對此進行改進。目前，該網站上沒有人有足夠的聲譽接受編輯。我很高興您建議通過評論進行編輯，或者編輯帖子，然後等到有人達到350代表。編輯完成後，我們可以刪除元註釋。
好！我編輯了答案；似乎畢竟是一種StackExchange-y的處理方式。
我也很喜歡Python，但是我覺得您不需要編寫腳本來完成如此瑣碎的事情-有許多命令行工具可以更快地完成此任務。
-1
@KirillG我通常會打開REPL，但我認為我仍然傾向於使用“ awk”（或類似的語言）。幾年前是一個不同的故事，但是自從改用Linux以來，學習如何使用“ awk”的基礎知識一直是我提高生產力最多的工作之一。但是可以肯定的是，我絕對同意這是否是更大的腳本或管道的一部分，我不敢稱之為`shell`。雖然我可能會使用`pysam`而非`biopython`，但這只是個人喜好:)
@KirillG我編寫了一個BioPython解決方案，該解決方案可以說更具可讀性/ Python風格，並且不會將整個輸入文件放入內存中https://bioinformatics.stackexchange.com/a/106/104
@Chris_Rands'的答案比我的要好得多，請投票！
謝謝Scott，儘管使用@KirillG's編輯您的答案也更像pythonic！但是，只有最新版本的BioPython允許將字符串作為`SeqIO.parse`的輸入（而不是文件句柄），在這種情況下，我認為BioPython從未明確關閉輸入文件句柄嗎？
嗯關於Biopython版本的要點-我僅在Python 2.7.11和3.5.1上測試了最新版本。我也對關閉文件句柄有些緊張，但是SeqIO.parse返回一個生成器，而不是一個文件，所以當我們到達StopIteration時，它不會自動關閉句柄嗎？ （免責聲明-我沒有任何證據！）
Chris_Rands
2017-05-18 13:04:25 UTC
view on stackexchange narkive permalink

這是 BioPython 的一種方法。 with 語句可確保關閉輸入和輸出文件句柄，並採用 lazy 方法，以便僅保存一條fasta記錄。一次存儲，而不是將整個文件讀入內存，這對於大型輸入文件是個壞主意。該解決方案不對序列ID長度或序列分佈的行數做任何假設：

 來自Bio import SeqIO，其中open（'sequences.fasta'）為in_f，open（' sequence.bed'，'w'）as out_f：記錄在SeqIO.parse（in_f，'fasta'）中：out_f.write（'{} \ t0 \ t {} \ n'.format（record.id，len （記錄）））
SmallChess
2017-05-19 10:14:21 UTC
view on stackexchange narkive permalink

我們有許多出色的答案！對於將來的用戶來說，這將是一個很好的參考。

我發現了自己在問的問題到底是什麼：

https：//www.biostars。 org / p / 191052 / 

  $ pip install pyfaidx $ faidx --transform bed test.fasta > test.bed

這是我要的單行命令。其他答案也可以，但是我想接受自己的答案。

pyfaidx的開發者在這裡。我正要寫下這個答案，並感嘆我沒有選擇一個更明顯或更令人難忘的軟件包名稱。
gringer
2017-05-18 06:04:49 UTC
view on stackexchange narkive permalink

如果FASTA序列全部在一行上，則以下perl一線應起作用：

  cat myseqs.fasta | perl -ne'if（/ ^ >（[^] +）/）{print $ 1} else {print“ 0”，length，“ \ n”}'

說明：

  • 如果行以'>'開頭，則將所有內容打印到第一個空格（但不要在行尾插入換行符）
  • 否則，打印“ 0”，然後是行的長度，然後是換行符


ⓘ
該問答將自動從英語翻譯而來。原始內容可在stackexchange上找到，我們感謝它分發的cc by-sa 3.0許可。
Loading...