Add parser to summarize diamond output

[olgabot]

For each query item, summarize the highest matching protein sequences as output by DIAMOND.

The query item differs depending on the input:

• if protein sequences provided --> each
• If bam + bed provided --> each interval
• If hashes provided --> each hash
• If differential hash expression (https://github.com/czbiohub/nf-predictorthologs/pull/26) --> each group + each hash differentially enriched in that group

Summarize the DIAMOND blast output, which looks like this:

 Fri 10 Apr - 09:58  ~/data_sm/tabula-microcebus/analyses/de-novo-orthology/kmermaid--muris-microcebus--v0
 olga@lrrr  head -n 3   mostly_marrow__unaligned/diamond/blastp/*.tsv | head -n 20
==> mostly_marrow__unaligned/diamond/blastp/hash-10140204645043044438__diamond__vertebrate_mammalian_concatenated.faa_db.tsv <==
A00111:82:H523JDMXX:1:1144:21766:19241  NP_573508.1     100.0   2.6e-07 58.5    NP_573508.1 alpha-hemoglobin-stabilizing protein [Mus musculus] 10090   Mus musculus    Eukaryota       Chordata
A00111:82:H523JDMXX:1:1144:21766:19241  XP_021024768.1  96.6    5.8e-07 57.4    XP_021024768.1 alpha-hemoglobin-stabilizing protein [Mus caroli]        10089   Mus caroli      Eukaryota       Chordata
A00111:82:H523JDMXX:1:1144:21766:19241  XP_021075064.1  89.7    8.3e-06 53.5    XP_021075064.1 alpha-hemoglobin-stabilizing protein [Mus pahari]        10093   Mus pahari      Eukaryota       Chordata

==> mostly_marrow__unaligned/diamond/blastp/hash-10347257497441110679__diamond__vertebrate_mammalian_concatenated.faa_db.tsv <==
A00111:82:H523JDMXX:1:2117:22064:10238  NP_573508.1     100.0   3.2e-08 61.6    NP_573508.1 alpha-hemoglobin-stabilizing protein [Mus musculus] 10090   Mus musculus    Eukaryota       Chordata
A00111:82:H523JDMXX:1:2117:22064:10238  XP_021075064.1  96.7    9.2e-08 60.1    XP_021075064.1 alpha-hemoglobin-stabilizing protein [Mus pahari]        10093   Mus pahari      Eukaryota       Chordata
A00111:82:H523JDMXX:1:2117:22064:10238  XP_021024768.1  96.7    1.2e-07 59.7    XP_021024768.1 alpha-hemoglobin-stabilizing protein [Mus caroli]        10089   Mus caroli      Eukaryota       Chordata

==> mostly_marrow__unaligned/diamond/blastp/hash-10513776695528344320__diamond__vertebrate_mammalian_concatenated.faa_db.tsv <==
A00111:82:H523JDMXX:1:1210:20021:21887  NP_573508.1     100.0   3.4e-07 58.2    NP_573508.1 alpha-hemoglobin-stabilizing protein [Mus musculus] 10090   Mus musculus    Eukaryota       Chordata

==> mostly_marrow__unaligned/diamond/blastp/hash-10858025809449178906__diamond__vertebrate_mammalian_concatenated.faa_db.tsv <==
A00111:82:H523JDMXX:1:1264:1759:7576    NP_573508.1     100.0   1.3e-09 66.2    NP_573508.1 alpha-hemoglobin-stabilizing protein [Mus musculus] 10090   Mus musculus    Eukaryota       Chordata
A00111:82:H523JDMXX:1:1264:1759:7576    XP_021024768.1  96.6    2.4e-08 62.0    XP_021024768.1 alpha-hemoglobin-stabilizing protein [Mus caroli]        10089   Mus caroli      Eukaryota       Chordata
A00111:82:H523JDMXX:1:1264:1759:7576    XP_021075064.1  96.6    1.2e-07 59.7    XP_021075064.1 alpha-hemoglobin-stabilizing protein [Mus pahari]        10093   Mus pahari      Eukaryota       Chordata

==> mostly_marrow__unaligned/diamond/blastp/hash-11316486562207166483__diamond__vertebrate_mammalian_concatenated.faa_db.tsv <==
A00111:82:H523JDMXX:1:1231:14543:17315  XP_021024768.1  100.0   6.8e-08 60.5    XP_021024768.1 alpha-hemoglobin-stabilizing protein [Mus caroli]        10089   Mus caroli      Eukaryota       Chordata

This summarizer could filter for NP-only sequences (https://github.com/czbiohub/nf-predictorthologs/issues/28)

[olgabot]

Here's some starter code:

filenames = glob.iglob(f'{outdir}/diamond/blastp/*.tsv')

dfs = []

DIAMOND_BLASTP_COLUMNS = ['read_id', 'subject_id', 'percent_identity', 'e_value', 'bitscore', 
           'subject_title', 'subject_taxid', 'subject_species', 'subject_kingdom', 
           'subject_superkingdom',
           'subject_phylum']

def read_diamond_blastp_output(filename):
    df = pd.read_csv(filename, sep='\t', names=DIAMOND_BLASTP_COLUMNS)
    return df

n = 0

for filename in tqdm(filenames):
    filesize = os.path.getsize(filename)
    if filesize > 0:
        n += 1
        basename = os.path.basename(filename)
        split = basename.split('__')
        utar_id = split[0]

        df = read_diamond_blastp_output(filename)
        df['region_id'] = utar_id
        dfs.append(df)

print(f'Number of regions with homology: {n}')
diamond_results = pd.concat(dfs)
print(diamond_results.shape)
diamond_results.head()

Extract gene predictions

Ignore the NP_01232 and species part of the description

pattern = '\d+.\d(.+)\[[\w ]+\]'
diamond_results.subject_title.head().str.extract(pattern)

Add more columns

Add columns for just the subject description, whether it's only a predicted gene or a "real" gene

diamond_results['subject_description_with_predicted'] = diamond_results.subject_title.str.extract(pattern)
diamond_results['is_predicted'] = diamond_results.subject_title.str.contains('PREDICTED')
diamond_results['subject_description'] = diamond_results.subject_description_with_predicted.str.split('PREDICTED: ').str[-1]
diamond_results.head()

If bed file is provided, compute most likely sequence using a weighted average

described here: https://github.com/czbiohub/nf-predictorthologs/issues/61, with most_likely_sequence function defined

diamond_results['e_value_inverse'] = 1/diamond_results['e_value']

diamond_predictions = diamond_results.groupby('region_id').apply(lambda x: most_likely_sequence(x, id_col='subject_description', weight_col='e_value_inverse'))
print(diamond_predictions.shape)
diamond_predictions.head()