consensus_insert_sequence.py

from Bio import SeqIO
from Bio import AlignIO
from Bio.Align import AlignInfo
from io import StringIO
from Bio.Align.Applications import MafftCommandline
from pathlib import Path
import argparse

# Functions to calculate a consensus sequence from SV inserted sequences
# alignment = AlignIO.read("results/insert_seqs/03-005-0130_03-5994_nl/ID_10002.aln", "clustal")

def calculate_percent_identity(alignment):
    alignment_length = alignment.get_alignment_length()
    total_identical_positions = 0 
    
    # Run through each position of the alignment
    for i in range(alignment_length):
        column = alignment[:, i]
        
        # Check if all sequences have the same character at this position
        if len(set(column)) == 1:
            total_identical_positions += 1 
    
    # Calculate and return the percent identity
    return (total_identical_positions / alignment_length) * 100

def get_consensus(msa_fasta):
    mafft_cline = MafftCommandline(input=msa_fasta, clustalout = True, maxiterate = 1000, localpair = True)
    stdout, stderr = mafft_cline()
    alignment = AlignIO.read(StringIO(stdout), "clustal")
    summary_align = AlignInfo.SummaryInfo(alignment)
    # Generate consensus with a consensus of 70% common residues at each position, and at least 2 sequences, otherwise, output N
    consensus = summary_align.dumb_consensus(threshold=0.7, ambiguous='N', require_multiple=2)

    # Output depth, alignment length, consenus and % identity
    perc_ident=calculate_percent_identity(alignment)

    return {'depth':len(alignment), 'length':alignment.get_alignment_length(), 'perc_identity':perc_ident, 'consensus':consensus.upper()}

def main():
    parser = argparse.ArgumentParser(description='Template Python3 Script')
    parser.add_argument('-i', '--input_ins_seqs', help = 'Inserted sequences generated by Savana, FASTA format', required = True)
    parser.add_argument('-n', '--insertion_id', help = 'Identifier of the inserion in the Savana VCF', required = True)
    parser.add_argument('-o', '--out_dir', default='.', help = 'Output directory')

    args = parser.parse_args()

    # Get consensus sequence
    consensus_dict=get_consensus(args.input_ins_seqs)

    # Write to file
    out_file = Path(args.out_dir, args.insertion_id + '.csv')
    with open(out_file, "w") as f:
        f.write(f'{args.insertion_id},{consensus_dict["depth"]},{consensus_dict["length"]},{consensus_dict["perc_identity"]},{consensus_dict["consensus"]}')

if __name__ == "__main__":
        main()