{"action":"created","issue":{"url":"https://api.github.com/repos/chapmanb/bcbio-nextgen/issues/859","labels_url":"https://api.github.com/repos/chapmanb/bcbio-nextgen/issues/859/labels{/name}","comments_url":"https://api.github.com/repos/chapmanb/bcbio-nextgen/issues/859/comments","events_url":"https://api.github.com/repos/chapmanb/bcbio-nextgen/issues/859/events","html_url":"https://github.com/chapmanb/bcbio-nextgen/issues/859","id":79444919,"number":859,"title":"RSEM is not using disambiguated reads but original inputs for transcriptome mapping","user":{"login":"schelhorn","id":5024959,"avatar_url":"https://avatars.githubusercontent.com/u/5024959?v=3","gravatar_id":"","url":"https://api.github.com/users/schelhorn","html_url":"https://github.com/schelhorn","followers_url":"https://api.github.com/users/schelhorn/followers","following_url":"https://api.github.com/users/schelhorn/following{/other_user}","gists_url":"https://api.github.com/users/schelhorn/gists{/gist_id}","starred_url":"https://api.github.com/users/schelhorn/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/schelhorn/subscriptions","organizations_url":"https://api.github.com/users/schelhorn/orgs","repos_url":"https://api.github.com/users/schelhorn/repos","events_url":"https://api.github.com/users/schelhorn/events{/privacy}","received_events_url":"https://api.github.com/users/schelhorn/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":4,"created_at":"2015-05-22T12:54:19Z","updated_at":"2015-06-23T07:38:46Z","closed_at":null,"body":"When running with `disambiguate` and `rsem`, is occured to me that RSEM is not using the disambiguated reads as inputs for the transcriptome mapping, but the original FASTQ files that still contain explant contamination. As a result, the RSEM results are not reliable in this `disambiguate` setting. This behavior is not transparent to the user (and also rather unexpected) and may therefore lead to false conclusions in downstream analyses.\r\n\r\nI suggest that either RSEM transcriptome mapping uses the disambiguated BAM file as an input (preferably directly without prior conversion to FASTQ, however the latter does not seem possible in the [default transcriptome mapper, bowtie2](https://github.com/chapmanb/bcbio-nextgen/blob/96542d8f5fa9ef3f19df78a45a2a12609bffa956/bcbio/pipeline/rnaseq.py#L70)), or that RSEM processing is disabled for the time being if disambiguation is set (by throwing an exception or skipping RSEM with a warning). \r\n\r\n@roryk, perhaps you could have a look since you seem to be the primary author of the RSEM functionality. Thank you."},"comment":{"url":"https://api.github.com/repos/chapmanb/bcbio-nextgen/issues/comments/114393448","html_url":"https://github.com/chapmanb/bcbio-nextgen/issues/859#issuecomment-114393448","issue_url":"https://api.github.com/repos/chapmanb/bcbio-nextgen/issues/859","id":114393448,"user":{"login":"schelhorn","id":5024959,"avatar_url":"https://avatars.githubusercontent.com/u/5024959?v=3","gravatar_id":"","url":"https://api.github.com/users/schelhorn","html_url":"https://github.com/schelhorn","followers_url":"https://api.github.com/users/schelhorn/followers","following_url":"https://api.github.com/users/schelhorn/following{/other_user}","gists_url":"https://api.github.com/users/schelhorn/gists{/gist_id}","starred_url":"https://api.github.com/users/schelhorn/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/schelhorn/subscriptions","organizations_url":"https://api.github.com/users/schelhorn/orgs","repos_url":"https://api.github.com/users/schelhorn/repos","events_url":"https://api.github.com/users/schelhorn/events{/privacy}","received_events_url":"https://api.github.com/users/schelhorn/received_events","type":"User","site_admin":false},"created_at":"2015-06-23T07:38:46Z","updated_at":"2015-06-23T07:38:46Z","body":"Alright, I came back to this issue and implemented what I believe should enable RSEM for disambiguated BAMs without breaking anything else. I would appreciate if @roryk or @chapmanb could have a look and integrate the code if appropriate. It seems that RSEM is not really used in bcbio anyway, but these fixes also re-enable the transcriptome BAM fir eXpress, so it should be generally useful.\r\n\r\n[bcbio.ngsalign.alignprep line 299](https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/ngsalign/alignprep.py#L299):\r\n\r\n```python\r\ndef _bgzip_from_bam(bam_file, dirs, config, is_retry=False, output_infix=''):\r\n``` \r\n[bcbio.ngsalign.alignprep line 310](https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/ngsalign/alignprep.py#L310):\r\n\r\n```python\r\n    out_file_1 = os.path.join(work_dir, \"%s%s-1.fq.gz\" % (os.path.splitext(os.path.basename(bam_file))[0]), output_infix)\r\n```\r\n[bcbio.pipeline.rnaseq line 56 ff](https://github.com/chapmanb/bcbio-nextgen/blob/master/bcbio/pipeline/rnaseq.py#L56) (The import belongs to the top of that file, of course):\r\n\r\n```python\r\nfrom bcbio.ngsalign import alignprep\r\ndef generate_transcript_counts(data):\r\n    \"\"\"Generate counts per transcript and per exon from an alignment\"\"\"\r\n    data[\"count_file\"] = featureCounts.count(data)\r\n\r\n    if dd.get_fusion_mode(data, False):\r\n        oncofuse_file = oncofuse.run(data)\r\n        if oncofuse_file:\r\n            data = dd.set_oncofuse_file(data, oncofuse_file)\r\n\r\n    if dd.get_rsem(data) and not dd.get_transcriptome_bam(data):\r\n        file1, file2 = None, None\r\n\r\n        if dd.get_disambiguate(data):\r\n            bam_path = data[\"work_bam\"]\r\n            fastq_paths = alignprep._bgzip_from_bam(bam_path, data[\"dirs\"], data[\"config\"], is_retry=False, output_infix='-transcriptome')\r\n            if len(fastq_paths) == 2:\r\n                file1, file2 = fastq_paths\r\n            else:\r\n                file1, file2 = fastq_paths[0], None\r\n        else:\r\n            file1, file2 = dd.get_input_sequence_files(data)\r\n\r\n        ref_file = dd.get_ref_file(data)\r\n        logger.info(\"RSEM was flagged to run, but the transcriptome BAM file \"\r\n                    \"was not found. Aligning to the transcriptome with bowtie2.\")\r\n        data = bowtie2.align_transcriptome(file1, file2, ref_file, data)\r\n    return [[data]]\r\n```"}}