Decoding Prolonged COVID-19: From Cluster Analysis to Pathway Mapping

bwa index Reference_genome.fasta

bwa mem Reference_genome.fasta sample_ID.fasta > sample_ID.sam

samtools view -@ 4 -Sb -o sample_ID.bam sample_ID.sam

freebayes -f reference.fasta input.bam > variants.vcf

bgzip filename.vcf

tabix filename.vcf

samtools faidx reference_seq.fasta MN908947.3 | bcftools consensus C1.vcf.gz > C1.fa

blastn -query C1.fa -subject reference_seq.fasta -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qseq sseq" -out C1_vs_refseq.out

import pandas as pd

def parse_blast_output(blast_output_file):
    # List to store SNP information
    snp_data = []

    with open(blast_output_file) as f:
        for line in f:
            # Skip header lines or empty lines
            if line.startswith("#") or not line.strip():
                continue

            cols = line.strip().split("\t")
            print(f"Number of columns: {len(cols)}")
            print(f"Columns: {cols}")

            # Check if the line has at least 13 columns (including sequences)
            if len(cols) < 13:  
                print(f"Skipping incomplete line: {line}")
                continue

            try:
                query_id = cols[0]
                subject_id = cols[1]
                percent_identity = float(cols[2])
                alignment_length = int(cols[3])
                mismatches = int(cols[4])
                gap_opens = int(cols[5])
                query_start = int(cols[6])
                query_end = int(cols[7])
                subject_start = int(cols[8])
                subject_end = int(cols[9])
                e_value = float(cols[10])
                bit_score = float(cols[11])
                query_seq = cols[12]
                subject_seq = cols[13]

                # Compare sequences and identify SNPs
                for i in range(min(len(query_seq), len(subject_seq))):
                    if query_seq[i] != subject_seq[i]:
                        query_pos = query_start + i
                        subject_pos = subject_start + i
                        snp = {
                            'Query ID': query_id,
                            'Subject ID': subject_id,
                            'Position': subject_pos,
                            'Reference Base': subject_seq[i],
                            'Variant Base': query_seq[i]
                        }
                        snp_data.append(snp)
                        print(f"SNP found at position {subject_pos}: {subject_seq[i]} -> {query_seq[i]}")
            except IndexError as e:
                print(f"Error processing line: {line}")
                print(f"Exception: {e}")
            except ValueError as e:
                print(f"Error converting to integer: {line}")
                print(f"Exception: {e}")

    # Convert the SNP data to a DataFrame
    df = pd.DataFrame(snp_data)
    # Create an output filename based on the input filename
    output_filename = blast_output_file.replace(".out", "_SNPs.xlsx")
    # Save the DataFrame to an Excel file
    df.to_excel(output_filename, index=False)

# Example usage
for i in range(1,57):
    parse_blast_output(f"C{i}_vs_ref.out")
    i += 1