在 GitHub 上編輯此頁面

序列清理器

描述

我想分享我使用 Biopython 清理序列的腳本。您應該知道分析不良數據會耗費 CPU 時間,而解釋不良數據的結果會耗費人力時間,因此進行預處理始終很重要。

讓我將我的腳本稱為 “Sequence_cleaner”,其主要概念是移除重複序列、移除過短的序列(使用者定義最短長度)以及移除含有過多未知核苷酸 (N) 的序列(使用者定義允許的 N 百分比),最後使用者可以選擇是否要將結果輸出到檔案或直接列印。

腳本

import sys
from Bio import SeqIO


def sequence_cleaner(fasta_file, min_length=0, por_n=100):
    # Create our hash table to add the sequences
    sequences = {}

    # Using the Biopython fasta parse we can read our fasta input
    for seq_record in SeqIO.parse(fasta_file, "fasta"):
        # Take the current sequence
        sequence = str(seq_record.seq).upper()
        # Check if the current sequence is according to the user parameters
        if (
            len(sequence) >= min_length
            and (float(sequence.count("N")) / float(len(sequence))) * 100 <= por_n
        ):
            # If the sequence passed in the test "is it clean?" and it isn't in the
            # hash table, the sequence and its id are going to be in the hash
            if sequence not in sequences:
                sequences[sequence] = seq_record.id
            # If it is already in the hash table, we're just gonna concatenate the ID
            # of the current sequence to another one that is already in the hash table
            else:
                sequences[sequence] += "_" + seq_record.id

    # Write the clean sequences

    # Create a file in the same directory where you ran this script
    with open("clear_" + fasta_file, "w+") as output_file:
        # Just read the hash table and write on the file as a fasta format
        for sequence in sequences:
            output_file.write(">" + sequences[sequence] + "\n" + sequence + "\n")

    print("CLEAN!!!\nPlease check clear_" + fasta_file)


userParameters = sys.argv[1:]

try:
    if len(userParameters) == 1:
        sequence_cleaner(userParameters[0])
    elif len(userParameters) == 2:
        sequence_cleaner(userParameters[0], float(userParameters[1]))
    elif len(userParameters) == 3:
        sequence_cleaner(
            userParameters[0], float(userParameters[1]), float(userParameters[2])
        )
    else:
        print("There is a problem!")
except:
    print("There is a problem!")

使用命令列,您應該執行

python sequence_cleaner.py input[1] min_length[2] min[3]

有 3 個基本參數

例如

python sequence_cleaner.py Aip_coral.fasta 10 10

供您參考:如果您不在意第二個和第三個參數,您只會移除重複的序列

問題、建議或改進

發送電子郵件至 genivaldo.gueiros@gmail.com