Sequence Post-processing Workflow for NNTC HIV env Sequences

public 1yr ago 0 bookmarks

View Workflow

nntc_hiv_brain_lymphoid_reservoirs — View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output, operation

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

Overview

Sequence post-processing of NNTC env sequences from brain and lymphoid tissues. SGA consensus sequences used as input to the workflow are found at resources/sequences.fa . Output from the geno2pheno coreceptor

Code Snippets

shell:
    "mafft {input} > {output}"

SnakeMake MAFFT API (EBI) From line 8 of rules/donor_align.smk

shell:
    "fasttree -nt {input} > {output}"

SnakeMake fasttree From line 8 of rules/fasttree.smk

script:
   "../scripts/functional_filter.jl"

SnakeMake From line 12 of rules/functional_filter.smk

script: "../scripts/hypermut2.py"

SnakeMake From line 5 of rules/hypermut.smk

run:
   from Bio import SeqIO, Seq, SeqRecord
   from collections import Counter
   import numpy as np

   # for simple consensus, ignoring ties
   def mode(collection):
      c = Counter(collection)
      return c.most_common(1)[0][0]

   recs = list(SeqIO.parse(str(input), "fasta"))
   seqs = [list(r.seq.upper()) for r in recs]
   seq_array = np.array(seqs)

   cons = []
   for i in range(seq_array.shape[1]):
      cons += list(mode(seq_array[:,i]))
   cons_record = SeqRecord.SeqRecord(Seq.Seq("".join(cons)), id = "CONSENSUS", description = "")
   SeqIO.write([cons_record] + recs, output[0], "fasta")

SnakeMake numpy Biopython From line 11 of rules/hypermut.smk

run:
   import pandas as pd
   dfs = [pd.read_table(f, sep = '\t') for f in input]
   pd.concat(dfs).to_csv(output[0], sep = '\t', index = False)

SnakeMake Pandas From line 35 of rules/hypermut.smk

shell:
    "mafft {input} > {output}"

SnakeMake MAFFT API (EBI) From line 8 of rules/mafft_align.smk

run:
    from Bio import SeqIO
    import sys
    import re
    import os

    out_recs = []
    out_dir = os.path.dirname(str(output))

    with open(str(input), 'r') as io:
        for r in SeqIO.parse(io, "fasta"):
            donor = re.match("NNTC_[0-9]{2}", r.id).group()
            if donor == wildcards.donor:
                out_recs.append(r)

    SeqIO.write(out_recs, str(output),"fasta")

SnakeMake Biopython From line 6 of rules/split.smk

import re
import sys
from scipy.stats import fisher_exact
import pandas as pd

inf = float('inf')

IUPAC_CODES = {
    'R': r'[RAG]',
    'Y': r'[YCT]',
    'B': r'[BYSKCGT]',
    'D': r'[DRWKAGT]',
    'H': r'[HYWMACT]',
    'V': r'[VRSMAGC]',
    'N': r'[NRYBDHVWSKMACGT]',
    'W': r'[WAT]',
    'S': r'[SCG]',
    'K': r'[KGT]',
    'M': r'[MAC]',
}

VALID_NA_PATTERN = re.compile(r'[NRYBDHVWSKMACGT]')


def expand_iupac(pattern):
    result = []
    for char in pattern:
        if char in IUPAC_CODES:
            result.append(IUPAC_CODES[char])
        else:
            result.append(char)
    return ''.join(result)


DEFAULT_PATTERNS = {
    'hypermut_from': re.compile(expand_iupac(r'^G(?=RD)')),
    'hypermut_to': re.compile(expand_iupac(r'^A(?=RD)')),
    'control_from': re.compile(expand_iupac(r'^G(?=YN|RC)')),
    'control_to': re.compile(expand_iupac(r'^A(?=YN|RC)')),
}


def fasta_reader(filename):
    with open(filename) as fp:
        header = None
        seq = []
        for line in fp:
            if line.startswith('#'):
                continue
            elif line.startswith('>'):
                if seq:
                    yield header, ''.join(seq).upper()
                header = line[1:].strip().split(' ')[0]
                seq = []
            else:
                seq.append(line.strip())
        if seq:
            yield header, ''.join(seq).upper()


def find_sites(seq, pattern, site_range):
    sites = []
    for offset in site_range:
        match = pattern.search(seq[offset:])
        if match:
            sites.append(offset)
    return sites


def get_comparable_sites(refseq, naseq):
    sites = []
    for offset, (ref, na) in enumerate(zip(refseq, naseq)):
        if not VALID_NA_PATTERN.match(ref):
            continue
        if not VALID_NA_PATTERN.match(na):
            continue
        sites.append(offset)
    return sites


def hypermut(refseq, naseq, patterns=DEFAULT_PATTERNS):
    comparable_sites = get_comparable_sites(refseq, naseq)
    potential_muts = find_sites(
        refseq, patterns['hypermut_from'], comparable_sites)
    potential_ctrls = find_sites(
        refseq, patterns['control_from'], comparable_sites)
    matched_muts = find_sites(
        naseq, patterns['hypermut_to'], potential_muts)
    matched_ctrls = find_sites(
        naseq, patterns['control_to'], potential_ctrls)
    num_potential_muts = len(potential_muts)
    num_matched_muts = len(matched_muts)
    num_potential_ctrls = len(potential_ctrls)
    num_matched_ctrls = len(matched_ctrls)
    try:
        oddsratio = (
            (num_matched_muts / num_potential_muts) /
            (num_matched_ctrls / num_potential_ctrls)
        )
    except ZeroDivisionError:
        oddsratio = inf
    _, p = fisher_exact([
        [num_matched_muts, num_potential_muts - num_matched_muts],
        [num_matched_ctrls, num_potential_ctrls - num_matched_ctrls]
    ], 'greater')
    return (
        num_matched_muts,
        num_potential_muts,
        num_matched_ctrls,
        num_potential_ctrls,
        oddsratio,
        p)


def main():
    fasta_filename = snakemake.input[0]
    sequences = list(fasta_reader(fasta_filename))
    _, refseq = sequences.pop(0)
    cols = ["ID", "n_matched", "n_potential", "n_matched_ctrl", "n_potential_cntrl", "oddsratio", "p"]
    res = []
    for naid, naseq in sequences:
        res.append([naid] + list(hypermut(refseq, naseq, DEFAULT_PATTERNS)))
    df = pd.DataFrame(res, columns = cols)
    df.to_csv(snakemake.output[0], sep='\t', index = False)


if __name__ == '__main__':
    main()