Missing identifiers allelematrix / callsets / samples

[matthijsbrouwer]

Allelematrix doesn't link to callsets

1484 from the 6737 unique callSetDbIds from the barley allelematrix endpoint seem not to be included in the callsets endpoint, for example:

 'SAMEA104650675',
 'SAMEA104634582',
 'SAMEA104630144',
 'SAMEA104648159',
 'SAMEA104650723',
 'SAMEA104629428',
  ...

5225 from the 5225 unique callSetDbIds from the wheat allelematrix endpoint seem not to be included in the callsets endpoint, for example:

 'SAMEA7824966',
 'SAMEA7825972',
 'SAMEA7828512',
 'SAMEA7826208',
 'SAMEA7821926',
 'SAMEA7827342',
  ...

Callsets / Samples

Callsets includes sampleDbId as an integer:

{
  "additionalInfo":{},
  "callSetDbId":"SAMEA14020287",
  "created":"2020-02-08",
  "externalReferences":[],
  "sampleDbId":3627152,
  "updated":"2020-02-08",
  "variantSetDbIds":[]
}

However, samples includes sampleDbId as string:

{
  "additionalInfo":null,
  "column":null,
  "externalReferences":[{"referenceId":"SAMEA14020287","referenceSource":""}],
  "germplasmDbId":"AB_07378",
  "observationUnitDbId":null,
  "plateDbId":null,
  "plateName":null,
  "programDbId":null,
  "row":null,
  "sampleBarcode":null,
  "sampleDbId":"3627152",
  "sampleDescription":"DB_ID: SAMEA14020287; sample volume: n/a; concentration: n/a; concentration: n/a; RNASeq-only: n/a; Specify number of lanes: 1; buffer: n/a; library protocol: GBS (narrow size range); Run type & workflow: n/a; Run type & workflow: Single read; flow cell type: SP; Unique dual Index: Meyer K,L,M; of index read cycles: 8 cycles; Average size (Tape): 412; size range: 400-600; size fractionation: BluePippin; REMARKS: n/a; REMARKS: n/a; spike-in: n/a; Sequenz: n/a; Read length: 100; Aussaat: n/a; library source: n/a; library strategy: n/a; library selection: n/a; sequencer/platform: Illumina NovaSeq 6000",
  "sampleGroupDbId":null,
  "sampleName":null,
  "samplePUI":"SAMEA14020287",
  "sampleTimestamp":"2020-02-08",
  "sampleType":"DNA",
  "studyDbId":null,
  "takenBy":null,
  "tissueType":"leaf",
  "trialDbId":null,
  "well":null}

Strictly speaking, this doesn't match.

Furthermore, possibly the data included under sampleDescription does fit under other BrAPI endpoints (?)

[Gouripriya5]

Callsets / Samples:

I have updated the sampleDbId in both samples and callsets to be in string format. They now match correctly.

[patrick-koenig]

@Gouripriya5
The missing callSetDbIds in the /callsets call are your or Stefans responsibility. Those missing BioSamples-IDs (=callsetDbIds) should be in the AGENT data warehouse. Can you investigate on this?

[Gouripriya5]

@matthijsbrouwer Thank you for bringing this issue to our attention.

I've reviewed the details you provided regarding the discrepancies between the AlleleMatrix and the CallSets endpoints. Upon reviewing the details, it appears there are two categories of AGENT sequence sets that we manage:

Core AGENT Sequences: These are consistent and all have bio sample IDs, ensuring that all sequence IDs, DbIDs, and CallSets IDs are available and correctly linked.

Additional Sequences from DivBrowse: The discrepancies you're seeing are due to the fact that these additional sequences were imported by Patrik from entire VCF files. They include sequences from sources such as GenBank 2.0. Genebank 2.0 do not have complete information, especially phenotypic data. They include germplasm info but lack phenotypic details.

The missing callSetDbIds you identified are likely from this second set, which is why they appear in the AlleleMatrix but not in our core CallSets endpoint.

Considering the absence of phenotypic information for these sequences, how do you suggest we handle these aspects moving forward?

Please let me know your thoughts on how we should proceed with integrating and managing these sequences. Thanks.

[matthijsbrouwer]

@patrick-koenig still some problems with the allematrix:

• selecting a specific callSetDbId or variantDbId or both does not seem to work
• also unknownString, sepPhased and sepUnphased don't seem to work
• please check all parameters, for example variantSetDbId

Also the datamatrix seems to be transposed, for example if I compare it with the output of the BrAPI-testserver:

• 2 variants and 3 callsets on BrAPI-testserver

[["0/0","1/0","1/0"],["0/0","0/0","0/0"]]

• 2 variants and 3 callsets on AGENT wheat implementation

[["-1/-1","-1/-1"],["0/0","0/0"],["-1/-1","-1/-1"]]

Also her -1 is used, but unknownString is the default "."

Furthermore, for the variants it would be nice if the chromosome can be included using at least the referenceName/DbId and possibly the reference genome using the referenceSetName/DbId

Finally, it would be nice if variantSet is used to distinguish barley/wheat, by using the commonCropName

[matthijsbrouwer]

@Gouripriya5 I understand your point, but it would be helpful to obtain this information directly from the BrAPI endpoints. At present, it can't be derived from either the allelematrix or callsets endpoints. We only know whether it is barley or wheat from the URL of the allelematrix endpoint, and even this will be obscured if it is behind the BrAPI-proxy solution. Therefore, incorporating these callSetDbIds into the callsets endpoint, linking them to at least some additional information, would be desirable. Without this, they might not be very useful.

[langeipk]

The arguments of @matthijsbrouwer makes perfectly sence. If there are genotyping data availale in form of sequences and VCF snp matrix, we shoul not hide them, althogh tehy are not reference-wise complete to all AGENT data domain, like phenotyping.
Thus, we asked our LIMS manager to include also sequence data that where sequenced in other projects, but where included in AGENT VCF file. The concrete criterion is: has a AGENT ID assigned and a BioSample ID was minted.
As Stephen informed me, AGENT IDs are assigned on dedicated wish by Nils to some selected Barley and Wheat acession that has a relevance for AGENT too. This list was submitted to Stephan bilateraly.

[Gouripriya5]

Allelematrix doesn't link to callsets

1484 from the 6737 unique callSetDbIds from the barley allelematrix endpoint seem not to be included in the callsets endpoint, for example:

 'SAMEA104650675',
 'SAMEA104634582',
 'SAMEA104630144',
 'SAMEA104648159',
 'SAMEA104650723',
 'SAMEA104629428',
  ...

5225 from the 5225 unique callSetDbIds from the wheat allelematrix endpoint seem not to be included in the callsets endpoint, for example:

 'SAMEA7824966',
 'SAMEA7825972',
 'SAMEA7828512',
 'SAMEA7826208',
 'SAMEA7821926',
 'SAMEA7827342',
  ...

Callsets / Samples

Callsets includes sampleDbId as an integer:

{
  "additionalInfo":{},
  "callSetDbId":"SAMEA14020287",
  "created":"2020-02-08",
  "externalReferences":[],
  "sampleDbId":3627152,
  "updated":"2020-02-08",
  "variantSetDbIds":[]
}

However, samples includes sampleDbId as string:

{
  "additionalInfo":null,
  "column":null,
  "externalReferences":[{"referenceId":"SAMEA14020287","referenceSource":""}],
  "germplasmDbId":"AB_07378",
  "observationUnitDbId":null,
  "plateDbId":null,
  "plateName":null,
  "programDbId":null,
  "row":null,
  "sampleBarcode":null,
  "sampleDbId":"3627152",
  "sampleDescription":"DB_ID: SAMEA14020287; sample volume: n/a; concentration: n/a; concentration: n/a; RNASeq-only: n/a; Specify number of lanes: 1; buffer: n/a; library protocol: GBS (narrow size range); Run type & workflow: n/a; Run type & workflow: Single read; flow cell type: SP; Unique dual Index: Meyer K,L,M; of index read cycles: 8 cycles; Average size (Tape): 412; size range: 400-600; size fractionation: BluePippin; REMARKS: n/a; REMARKS: n/a; spike-in: n/a; Sequenz: n/a; Read length: 100; Aussaat: n/a; library source: n/a; library strategy: n/a; library selection: n/a; sequencer/platform: Illumina NovaSeq 6000",
  "sampleGroupDbId":null,
  "sampleName":null,
  "samplePUI":"SAMEA14020287",
  "sampleTimestamp":"2020-02-08",
  "sampleType":"DNA",
  "studyDbId":null,
  "takenBy":null,
  "tissueType":"leaf",
  "trialDbId":null,
  "well":null}

Strictly speaking, this doesn't match.

Furthermore, possibly the data included under sampleDescription does fit under other BrAPI endpoints (?)

Hello Matthijs,

I wanted to follow up on the issue#7 you raised regarding the 5225 unique CallSetDbIds from the wheat allele matrix and 1484 from 6737 from the barley endpoint not being included in the callsets endpoint. You mentioned several examples, such as:

SAMEA7824966
SAMEA7825972
SAMEA7828512
SAMEA7826208
SAMEA7821926
SAMEA7827342

Based on your feedback, I conducted a search using the following URL:
https://divbrowse.ipk-gatersleben.de/agent_wheat/brapi/v2/allelematrix?dimensionCallSetPageSize=5225

However, in the results, I was unable to find some of the CallSetDbIds you listed, such as SAMEA7825972. When I expanded the search using a larger page size (100,000) with the URL:
https://divbrowse.ipk-gatersleben.de/agent_wheat/brapi/v2/allelematrix?dimensionCallSetPageSize=100000 , I was able to locate those CallSetDbIds and also same for barley.

Could you please clarify how you identified the 5225 unique CallSetDbIds in wheat and the 1484 from 6737 unique CallSetDbIds in barley? It would be helpful to understand your method so I can investigate the discrepancy further and ensure the endpoints are consistent.

[Gouripriya5]

Hi Priya,

Attached a jupyter notebook with some tests to (re)construct these problems

Best,
Matthijs
20240724_test.json

[matthijsbrouwer]

below some code to test

import requests

url_barley = "https://divbrowse.ipk-gatersleben.de/agent_barley/brapi/v2"
url_wheat = "https://divbrowse.ipk-gatersleben.de/agent_wheat/brapi/v2"
url_all = "https://agent.ipk-gatersleben.de/genotyping/brapi/v2"

def brapiRequest(baseUrl,call,params={}):
    try:
        headers = {"Accept": "application/json"}
        url = "{}/{}".format(baseUrl,call)
        fullUrl = "{}?{}".format(url,"&".join(["{}={}".format(x,y) for x,y in params.items()]))
        response = requests.get(url, params=params, headers=headers)
        try:
            if response.ok:
                return response.json(), response.status_code, None, fullUrl
            else:
                return None, response.status_code, response.text, fullUrl
        except:
            return None, 500, response.text, fullUrl
    except Exception as e:
        return None, 500, "error: {}".format(str(e)), None

def getAllelematrixVariantsCallsets(url):
    variants = []
    callsets = []
    totalVariantsPages = 1
    totalCallsetsPages = 1
    variantsPage = 0
    callsetsPage = 0
    while variantsPage<totalVariantsPages:
        (response, responseCode, responseError, responseUrl) = brapiRequest(url, "allelematrix", 
                                       {"preview":"true", 
                                        "dimensionVariantPage": variantsPage,
                                        "dimensionVariantPageSize": 10000,
                                        "dimensionCallSetPageSize": 1})
        for item in response["result"]["pagination"]:
            if item["dimension"]=="VARIANTS":
                totalVariantsPages = item["totalPages"]
        variants = variants + response["result"]["variantDbIds"]
        # print("variants",variantsPage,totalVariantsPages)
        variantsPage+=1
    while callsetsPage<totalCallsetsPages:
        (response, responseCode, responseError, responseUrl) = brapiRequest(url, "allelematrix", 
                                       {"preview":"true", 
                                        "dimensionCallSetPage": callsetsPage,
                                        "dimensionCallSetPageSize": 10000,
                                        "dimensionVariantPageSize": 1})
        for item in response["result"]["pagination"]:
            if item["dimension"]=="CALLSETS":
                totalCallsetsPages = item["totalPages"]
        callsets = callsets + response["result"]["callSetDbIds"]
        # print("callsets",callsetsPage,totalCallsetsPages)
        callsetsPage+=1
    return variants,callsets

allelematrix_variants_barley, allelematrix_callsets_barley = getAllelematrixVariantsCallsets(url_barley)
allelematrix_variants_wheat, allelematrix_callsets_wheat = getAllelematrixVariantsCallsets(url_wheat)

def getAllIdentifiers(url,call,identifier):
    page = 0
    totalPages = 1
    results = []
    while page<totalPages:
        (response, responseCode, responseError, responseUrl) = brapiRequest(url, call, 
                                       {"page": page, "pageSize": 10000})
        for item in response["result"]["data"]:
            results.append(item[identifier])
        totalPages = response["metadata"]["pagination"]["totalPages"]
        page+=1
    return results

variants_barley = getAllIdentifiers(url_barley,"variants","variantDbId")
variants_wheat = getAllIdentifiers(url_wheat,"variants","variantDbId")
#expect empty sets
print(set(allelematrix_variants_barley).difference(variants_barley))
print(set(allelematrix_variants_wheat).difference(variants_wheat))
print(set(variants_barley).difference(allelematrix_variants_barley))
print(set(variants_wheat).difference(allelematrix_variants_wheat))

callsets = getAllIdentifiers(url_all,"callsets","callSetDbId")

print("wheat: allelematrix %s callsets, %s are in callsets, %s not in callsets" % (
    len(set(allelematrix_callsets_wheat)),
    len(set(allelematrix_callsets_wheat).intersection(callsets)),
    len(set(allelematrix_callsets_wheat).difference(callsets))
))
print("barley: allelematrix %s callsets, %s are in callsets, %s not in callsets" % (
    len(set(allelematrix_callsets_barley)),
    len(set(allelematrix_callsets_barley).intersection(callsets)),
    len(set(allelematrix_callsets_barley).difference(callsets))
))
print("%s callsets not in allelematrix wheat or barley" % (
    len(set(callsets).difference(allelematrix_callsets_wheat).difference(allelematrix_callsets_barley))
))
print(sorted(set(allelematrix_callsets_wheat).difference(callsets))[0:10])
print(sorted(set(allelematrix_callsets_barley).difference(callsets))[0:10])
print(sorted(set(callsets).difference(allelematrix_callsets_wheat).difference(allelematrix_callsets_barley))[0:10])

#check examples from e-mail
entries = ["SAMEA7824966", "SAMEA7825972", "SAMEA7828512", "SAMEA7826208", "SAMEA7821926", "SAMEA7827342"]
for item in entries:
    print(item in allelematrix_callsets_wheat, item in callsets)

[Gouripriya5]

@matthijsbrouwer Thanks for the code.

Issue Summary:
Previously, there was a discrepancy where 1,484 out of 6,737 unique callSetDbIds from the barley allelematrix endpoint and all 5,225 unique callSetDbIds from the wheat allelematrix endpoint were not included in the callsets endpoint.

Resolution:
This issue has now been fully resolved. All callSetDbIds are correctly included across the relevant endpoints. The data is now synchronized, and you can verify that the callSetDbIds are present in both the allelematrix and callsets endpoints.

If you encounter any further discrepancies, please feel free to inform.

[langeipk]

referenceSource is empty: "referenceSource": ""
But because we deliver "referenceId" it should be filled. @Gouripriya5 Please check the SQL code and add a correct contant here to represent BioSample database as reference in accordance to BrAPI format for this attribute

[langeipk]

studydb id results in Stringt but should be the coreespanding numeric IDs frm study endpoint. Now we see e.g. "studyDbId": "Bridge_project" Please make sure that all projects that are devilered by sample endpoint are delivered from study endpoint