views.py 12.6 KB
Newer Older
1
# # Copyright 2013-2017 Mathieu Courcelles
2
# # Mike Tyers's lab / IRIC / Universite de Montreal 
Mathieu Courcelles's avatar
Mathieu Courcelles committed
3 4 5


# Import standard libraries
6
from math import floor
7
import numpy as np
8
import os
9
from scipy.stats import mannwhitneyu
Mathieu Courcelles's avatar
Mathieu Courcelles committed
10 11

# Import Django related libraries
12
from django.conf import settings
Mathieu Courcelles's avatar
Mathieu Courcelles committed
13
from django.contrib.auth.decorators import login_required
14
from django.http import Http404, HttpResponse
15
from django.shortcuts import render
16
from django.utils.safestring import mark_safe
Mathieu Courcelles's avatar
Mathieu Courcelles committed
17

18 19
# Import BioPython libraries
from Bio.PDB.PDBParser import PDBParser
Mathieu Courcelles's avatar
Mathieu Courcelles committed
20 21

# Import project libraries
22 23 24 25 26 27
from .parser.exception import InvalidFileFormatException
from .queryset_operation import clpeptide_set_2_protein_set
from .models import CLPeptide
from .pdb_structure import (compute_cl_distance,
                            randomDistance,
                            )
28
from numpy.lib.function_base import average
29

30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
# Dictionary of peptide modifications with associated symbol
modification_name_symbol = {
    # 'Carbamidomethyl': 'cm',
    # 'Cysteinyl': 'cy',
    # 'Deamidated': 'd',
    # 'Dimethyl': 'dme',
    # 'Methyl': 'me',
    # 'Oxidation': 'ox',
    # 'Phospho': 'ph',
    # 'TMT6plex': 'tmt'
    '[15.99]': 'ox',
    '[57.02]': 'cm',
    '[156.08]': 'dss',
}


Mathieu Courcelles's avatar
Mathieu Courcelles committed
46 47 48 49 50 51 52 53 54 55
def alive(request):
    """
    Returns HttpResponse with status 204
    :param request: HttpRequest
    :return: HttpResponse
    """

    return HttpResponse(status=204)


56 57 58 59 60 61 62 63 64 65 66 67 68
def format_peptide_seq_with_mod(peptide):
    """
    Formats peptide sequence string with modification for spectrum viewer.
    :param peptide: Peptide object
    :return: Peptide sequence string with modification symbols.
    """

    for modification, symbol in modification_name_symbol.items():
        peptide = peptide.replace(modification, symbol)

    return peptide


69 70 71 72 73 74 75 76 77 78 79 80
@login_required
def clpeptide_msms(request, pk):
    """
    This view displays PSM MS/MS spectrum.
    :param request: HttpRequest object
    :param pk: Peptide primary key integer
    :return: HttpResponse
    """

    clpeptide = CLPeptide.objects.get(pk=pk)

    # Get MGF file
81 82 83
    try:
        run_name = clpeptide.run_name + '.mgf'

84
        peak_list_object = PeakList.objects.filter(file__endswith=run_name).last()
85 86 87 88 89 90

    except PeakList.DoesNotExist:
        raise Http404('No Peak list entry for %s.' % run_name)

    msms_file_name = os.path.join(settings.MEDIA_ROOT, str(peak_list_object.file))
    parser = MGFParser(msms_file_name)
91 92

    # Get peak list
93 94
    try:
        scan = parser.get_scan_at_index(clpeptide.scan_file_index)
95

96 97
    except IOError:
        raise Http404('Spectrum file not found (%s)' % msms_file_name)
98 99 100 101 102 103 104 105 106 107 108 109 110 111

    for i, value in enumerate(scan['peak_list']):
        scan['peak_list'][i] = '%s %s' % (value[0], value[1])

    peak_list = '\n'.join(scan['peak_list'])

    validation_disabled = ''

    if clpeptide.validated == '':
        clpeptide.validated = None

    context = {
        'myTolerance': 10,
        'myPrecursorZ': clpeptide.precursor_charge,
112 113
        'myPeptide': format_peptide_seq_with_mod(clpeptide.peptide1),
        'myPeptide2': format_peptide_seq_with_mod(clpeptide.peptide2),
114 115
        'myPeaklist': peak_list,
        'clpeptide': clpeptide,
116 117
        'validation_disabled': validation_disabled,
        'linked_formula': clpeptide.dataset.first().cross_linker.linked_formula
118 119 120 121 122 123
    }

    return render(request,
                  'spectrum_viewer.html',
                  context)

124

125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
@login_required
def clpeptide_validation(request, pk, value):
    """
    This view updates validation status of a peptide spectrum match.
    :param request: HttpRequest object
    :param pk: Peptide primary key integer
    :param value: Validation character
    :return: HttpResponse object
    """

    clpeptide = CLPeptide.objects.get(pk=pk)

    json_data = '{"updated": 0}'

    if value in ['G', 'Y', 'R']:
        clpeptide.validated = value
        clpeptide.save()

        json_data = '{"updated": 1}'

    return HttpResponse(json_data, content_type='application/json')
Mathieu Courcelles's avatar
Mathieu Courcelles committed
146

147

Mathieu Courcelles's avatar
Mathieu Courcelles committed
148
@login_required
149
def xiNET_view(request, clpeptide_set, form):
Mathieu Courcelles's avatar
Mathieu Courcelles committed
150 151 152 153
    """
    This function generates a view to display cross-links with xiNET.
    """

154
    clpeptide_set = clpeptide_set.filter(not_decoy=True, cross_link=True)
155

156 157
    quant_set = QuantificationFC.objects.filter(quantification=form.cleaned_data['quantification'])

158 159
    # Prepare protein set
    protein_set = clpeptide_set_2_protein_set(clpeptide_set)
160

161 162 163
    # Find absolute max if none is specified
    limit = 0
    if form.cleaned_data['fold_change_limit'] is None:
164

165
        for quant in quant_set:
166

167
            absoluate_value = abs(quant.fold_change)
168

169 170 171 172
            if absoluate_value > limit:
                limit = absoluate_value
    else:
        limit = form.cleaned_data['fold_change_limit']
173

174 175 176 177 178 179 180 181
    c = {'clpeptide_set': clpeptide_set,
         'color_scheme': form.cleaned_data['color_scheme'],
         'fold_change_limit': limit,
         'protein_set': protein_set,
         'quant_set': quant_set,
         'hide_not_quantified': form.cleaned_data['hide_not_quantified'],
         'min_score': form.cleaned_data['min_score'],
         'max_score': form.cleaned_data['max_score']}
182

183
    return render(request, 'xiNET.html', c)
184

185 186

@login_required
187
def jsmol_view(request, clpeptide_set, form=None):
188 189 190
    """
    This function generates a view to display cross-links with JSMol.
    """
191

192
    pdb = form.get_PDB()
193

194 195
    ukey = dict()
    ukey['proteinPpos'] = True
196 197 198

    try:
        clpeptide_set, alignments = compute_cl_distance(clpeptide_set,
199 200 201 202
                                                        pdb,
                                                        form.cleaned_data['protein_identity'],
                                                        form.cleaned_data['peptide_identity'],
                                                        unique_key=ukey)
203 204
    except InvalidFileFormatException as e:
        return HttpResponse('<h1>%s</h1>' % e)
205

206
    cross_links = []
207

208 209
    cl_count = len(clpeptide_set)
    mapped_cl_count = 0
210

Mathieu Courcelles's avatar
Mathieu Courcelles committed
211
    distance_list = []
212
    min_distance_list = []
213

214 215
    # Quant
    quant_set = QuantificationFC.objects.filter(quantification=form.cleaned_data['quantification'])
216

217 218 219
    # Find absolute max if none is specified
    limit = 0
    quant_clpep_key = {}
220

221
    for quant in quant_set:
222

223
        absoluate_value = abs(quant.fold_change)
224

225 226
        if absoluate_value > limit:
            limit = absoluate_value
227

228 229 230 231
        quant_clpep_key[quant.clpeptide.pk] = True

    if form.cleaned_data['fold_change_limit'] is not None:
        limit = form.cleaned_data['fold_change_limit']
232

233
    clink_residues_dict = dict()
234

235
    for clpep in clpeptide_set:
236

237 238 239
        if form.cleaned_data['hide_not_quantified'] and clpep.pk not in quant_clpep_key:
            continue

240 241
        if clpep.fs_prot1_id is None:
            continue
242

243
        protein_1 = '%s [%s]' % (clpep.fs_prot1_id.identifier, clpep.peptide_position1 + clpep.pep1_link_pos - 1)
244

245 246 247
        protein_2 = ''
        if clpep.fs_prot2_id is not None:
            protein_2 = '%s [%s]' % (clpep.fs_prot2_id.identifier, clpep.peptide_position2 + clpep.pep2_link_pos - 1)
248

249
        j = False;
250

251
        min_distance = 1000000000000
252

253
        for cldistance in clpep.distances:
254

255 256
            if cldistance.distance != '':
                j = True
Mathieu Courcelles's avatar
Mathieu Courcelles committed
257
                distance_list.append(float(cldistance.distance))
258 259
                clink_residues_dict[cldistance.residue_1_str] = cldistance.residue_1
                clink_residues_dict[cldistance.residue_2_str] = cldistance.residue_2
260

261 262
                if min_distance > float(cldistance.distance):
                    min_distance = float(cldistance.distance)
263 264 265

            cross_links.append(
                '{pk: "%s", prot1: "%s", prot2: "%s", res1: "%s", res2: "%s", score: "%.2f", distance: "%s", ident1: "%s", ident2: "%s"}' % (
266 267 268 269 270 271 272 273 274
                    clpep.pk,
                    str(protein_1),
                    str(protein_2),
                    cldistance.residue_1_str,
                    cldistance.residue_2_str,
                    clpep.match_score,
                    cldistance.distance,
                    cldistance.identity_1,
                    cldistance.identity_2,
275
                ))
276 277 278

        if min_distance != 1000000000000:
            min_distance_list.append(min_distance)
279

280 281
        if j:
            mapped_cl_count += 1
282 283 284 285 286 287

    pdb.clink_residues_list = clink_residues_dict.values()

    # Check distance distribution of cross-links to see if it is random
    p_values = []
    draw_count = 1000
288

289 290
    p_val_threshold = 0.05
    lower_count = 0
291

292
    s_mean = np.mean(min_distance_list)
293

294
    for x in range(0, draw_count):
295

296 297
        random_distance_list = randomDistance(pdb, form.cleaned_data['cross_linked_residues'],
                                              len(min_distance_list))
298

Mathieu Courcelles's avatar
Mathieu Courcelles committed
299
        if len(min_distance_list) > 0 and len(random_distance_list) > 0:
300

Mathieu Courcelles's avatar
Mathieu Courcelles committed
301 302 303 304 305
            try:
                p_values.append(mannwhitneyu(min_distance_list, random_distance_list)[1] * 2)
                r_mean = np.mean(random_distance_list)
                if p_values[-1] <= p_val_threshold:
                    lower_count += 1
306

Mathieu Courcelles's avatar
Mathieu Courcelles committed
307 308 309
            except ValueError:
                p_values.append(1)
        else:
310
            p_values.append(1)
311

312 313 314
    min_p_value = min(p_values)
    max_p_value = max(p_values)
    avg_p_value = average(p_values)
315

316
    p_value_str = 'p-value <= %.1e' % max_p_value
317 318
    ratio = lower_count / float(draw_count)

319 320
    if ratio < 1:
        p_value_str = 'p-value <= 0.05 for %d %% of the draws.' % (ratio * 100)
321 322

    # Get compound details
323 324 325
    parser = PDBParser()
    structure = parser.get_structure('tmp', pdb.file.file.name)
    pdb.compound = structure.header['compound']
326

327 328 329
    chain_list = list()
    for chain in structure.get_chains():
        chain_list.append(chain.id)
330

331
    molecules = list()
332

333 334 335
    for key in pdb.compound:
        molecules.append([key, pdb.compound[key].get('molecule', ''), pdb.compound[key].get('chain', '')])

Mathieu Courcelles's avatar
Mathieu Courcelles committed
336
    # Histogram
337
    all_bin_values, all_bin_labels = ([], [])
338 339
    min_bin_values, min_bin_labels = ([], [])
    rand_bin_values, rand_bin_labels = ([], [])
340

341
    floor_all_distance_list = [0, ]
342

343
    if len(distance_list) != 0:
344
        floor_all_distance_list = [floor(x) for x in distance_list]
345
        all_bin_values, all_bin_labels = np.histogram(floor_all_distance_list, bins=np.arange(0,
346 347 348
                                                                                              max(
                                                                                                  floor_all_distance_list) + 4,
                                                                                              2))
349 350
    if len(min_distance_list) != 0:
        floor_min_distance_list = [floor(x) for x in min_distance_list]
351
        min_bin_values, min_bin_labels = np.histogram(floor_min_distance_list, bins=np.arange(0,
352 353 354 355
                                                                                              max(
                                                                                                  floor_all_distance_list) + 4,
                                                                                              2))

356 357 358
    if len(random_distance_list) != 0:
        floor_rand_distance_list = [floor(x) for x in random_distance_list]
        rand_bin_values, rand_bin_labels = np.histogram(floor_rand_distance_list, bins=np.arange(0,
359 360 361
                                                                                                 max(
                                                                                                     floor_all_distance_list) + 4,
                                                                                                 2))
362

363 364 365 366 367
    # Replace fasta_id with identifier
    for key, value in alignments.items():
        seq = FastaDb_Sequence.objects.filter(pk=key)
        identifier = seq[0].identifier
        alignments[identifier] = alignments.pop(key)
368

369
    # Ready
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
    c = {'cross_links_array': mark_safe("[" + ','.join(cross_links) + "]"),
         'pdb': pdb,
         'cl_count': cl_count,
         'mapped_cl_count': mapped_cl_count,
         'viewer': form.cleaned_data['viewer'],
         'molecules': molecules,
         'chain_list': chain_list,
         'alignments': alignments,
         'all_bin_values': all_bin_values,
         'min_bin_values': min_bin_values,
         'rand_bin_values': rand_bin_values,
         'bin_labels': all_bin_labels,
         'color_scheme': form.cleaned_data['color_scheme'],
         'fold_change_limit': limit,
         'quant_set': quant_set,
         'min_val': 0,
         'draw_count': draw_count,
         'p_value': p_value_str,
         'min_score': form.cleaned_data['min_score'],
         'max_score': form.cleaned_data['max_score']
         }

    return render(request, 'jsmol.html', c)