環比IDS兩個FASTA文件

我有多個序列環比IDS兩個FASTA文件

cat file1.fasta 
>1 
ACGTCGAT 
>2 
ACTTTATT 
>3 
ACGGGG 

cat file2.fasta 
>1 
CCGGAGC 
>2 
TGTCAGTC 
>3 
CTACGTCTT

2個FASTA文件我也有標識每個FASTA文件，我想使用的ID來提取特定序列的列表，使2序列fasta，然後執行一些操作（對齊，計算距離）。

列表：

cat file1.list 
1 
3 
cat file2.list 
2 
1

在現實中這些FASTA文件，並列出成千上萬的長

我想遍歷在列表中的每一行以提取匹配的FASTA文件序列/線特定的id /行，然後將每個文件中的fasta序列組合到可以對齊的兩個序列fasta文件中。基本上，我希望每個fasta序列與它的「pair」成對對齊。

因此，根據此處的示例和列表ID順序，我想將file1.fasta中的fasta序列1與file2.fasta中的fasta序列2配對，然後移至下一對（來自file1的序列3）。 fasta，以及來自file2.fasta的序列1等）。基於id提取fasta序列相對比較簡單（有幾種方法），但其中一個是faOneRecord，它將輸入的fasta文件作爲要提取的fasta文件，然後輸入要查找的記錄/ id，然後返回fasta序列和頭：

faOneRecord <in.fa> <recordName>

因此，第一循環結束後，我會根據ID列表上創建這個文件：

>1 
ACGTCGAT 
>2 
TGTCAGTC

等。

我認爲這樣做相對容易，但我似乎無法達到那裏。然後，一旦我製作了兩個fasta序列，每個循環，我想對齊並獲得距離估計值，打印到文件並轉到下一個循環。其餘部分可能需要一些工作，並需要特定的程序，但我需要幫助，只需生成2個序列fasta提取/循環通過ID。

我想主要的問題是如何循環的ID，然後管那些ID作爲參數傳遞到faOneRecord命令

這可能是太具體了，如果讓我道歉，但對如何獲得任何想法開始將是有益的，非常感謝。

來源

2017-02-20 user95146

你真的需要使用bash嗎？這是一個可以用一種語言解決的問題，在這種語言中，您可以打開兩個fasta文件並將序列保存到一個數組，以及可以並行迭代列表文件的位置。在bash中這兩者都非常棘手。 – Marian

它不需要在bash中......我只需要最終在bash中調用程序。 python（biopython）解決方案可能更直接，但我不太熟練。你有建議如何在python中啓動它？ – user95146

下面是一個python解決方案的（不完整）草圖。正如我在評論中所說的，有兩個步驟：

首先，讀取數組中的兩個文件。如果你確定他們是完全按照你的榜樣，你可以忽略>x線：

fasta1 = [''] # make sure the first item is saved to fasta1[1], not fasta[0] 
for line in open('file1.fasta'): 
    if not line.startswith('>'): 
     fasta1.append(line.strip())

的for line in open()剛剛打開文件，並遍歷其行。

對file2執行相同操作。然後，你可以交替閱讀list文件，拿到號碼的開出並打印匹配序列：

for l1, l2 in zip(open('file1.list'), open('file2.list')): 
    print(fasta1[int(l1)]) 
    print(fasta2[int(l1)])

zip需要兩個文件，並在並行讀取它們中，從而第一時間執行循環，l1和l2分別包含第一行file1.list和file2.list;第二次，它是每個的第二行，等等。

來源

2017-02-20 01:11:13 Marian

謝謝。這讓我開始了。一個問題 - fasta序列的包裝...實際上，這些fasta序列中的很多序列長度爲幾百個到幾千bp，並且序列包裝在一行中。我用包裹的fasta測試了腳本，它只打印第一行...有任何方法來處理這個問題？也許biopython解析器是需要的？或者我可以線性化序列.... – user95146

如果biopython有解析這些文件的函數，請使用它。我提出的代碼只是一個快速入侵。爲了完成，我將如何處理純Python中的行：將行（'.strip（）'ed）存儲在第二個列表中，當遇到'> n'行時，'''.join（）'所有這些行並將它們追加到'fasta1'。 – Marian

謝謝。如果在'list'文件中缺少ID，我們怎麼能錯誤處理'list Index outof Range'？我需要錯誤地處理每個索引調用嗎？例如'
for l1，l2 in zip（open（'file1。清單 '），開放（' file2.list'））：嘗試： A = fasta1 [INT（L1）除了IndexError，E：打印「GGG」' continue' – user95146

環比IDS兩個FASTA文件

回答

相關問題