1
$\begingroup$

I am trying snakemake for my mRNA-seq pipeline and I got stuck in a simple thing.

from os.path import join


RSEM_INDEX = "/app/ref/ensembl/human/idx/GRCh37"
SAMPLES, = glob_wildcards("Align/{sample}/Aligned.toTranscriptome.out.bam")

print("---Sample Names---")

for i in SAMPLES:
    print(i)

rule all:
    input:
        expand("Expression/{sample}/rsem.genes.results", sample=SAMPLES)


rule expression:
    input:
        bam = join("Align", "{wildcards.sample}}", "Aligned.toTranscriptome.out.bam")
    params:
        rsem_index = RSEM_INDEX
    output:
        join("Expression", "{wildcards.sample}", "rsem.genes.results")
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end Align/{sample}/Aligned.toTranscriptome.out.bam {params.rsem_index} Expression/{wildcards.sample}/rsem.genes.results >& Expression/{wildcards.sample}/rsem.log
        """

Here is an error I got

    MissingInputException in line 12 of /app/src/mRNAseq/TBD170049_20170320_RNAseq/quantify.snakemake:
Missing input files for rule all:
Expression/013_clamp_0m_buffy/rsem.genes.results
Expression/12_27_11_clamp_buffy_200/rsem.genes.results
Expression/010_200_buffy/rsem.genes.results
Expression/010_0_buffy/rsem.genes.results
Expression/12_27_11_clamp-p_buffy_0/rsem.genes.results
Expression/013_clamp_200m_buffy/rsem.genes.results
Expression/11_15_006_0m_buffy/rsem.genes.results
Expression/11_15_006_200m_buffy/rsem.genes.results

So, I thought rule all input for your final files you want to make, but it's looking for the final files that I want to make.

I would like to get some help. would you help me please? Thank you in advance!

My file trees.

drwxrwxr-x. 10 corea corea  4096 Sep 26 20:54 Align
-rw-rw-r--.  1 corea corea  1295 Sep 27 10:19 align.snakemake
-rw-rw-r--.  1 corea corea  1823 Sep 27 09:38 align.snakemake.test
-rw-rw-r--.  1 corea corea 31055 Sep 26 08:55 out
-rw-rw-r--.  1 corea corea   869 Sep 27 12:34 quantify.snakemake
drwxrwxr-x.  2 corea corea  4096 Sep 22 20:50 Rawdata
(mRNAseq) [corea@Biocore76 TBD170049_20170320_RNAseq]$ tree Align/
Align/
├── 010_0_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 010_200_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 013_clamp_0m_buffy
$\endgroup$
  • 2
    $\begingroup$ cross posted: biostars.org/p/274858 and stackoverflow.com/questions/46439175 $\endgroup$ – Pierre Sep 27 '17 at 6:29
  • $\begingroup$ I am sorry. I will keep that rule in mind and it will never happen! $\endgroup$ – jk_kim Sep 27 '17 at 23:14
  • $\begingroup$ @jk_kim I don't think cross-posting is wrong per se, provided this is clearly acknowledged and that links are given so that readers of one site know there may be interesting information on the other site. $\endgroup$ – bli Oct 27 '17 at 13:10
3
$\begingroup$

The following (or something close) should work:

rule expression:
    input:
        "Align/{sample}/Aligned.toTranscriptome.out.bam"
    params:
        rsem_index = RSEM_INDEX
    output:
        "Expression/{sample}/rsem.genes.results"
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end {input} {params.rsem_index} {output} >& Expression/{wildcards.sample}/rsem.log
        """

BTW, the original issue was probably due to the }} in your input line.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the comment! no more cross posting! I will keep that rule in mind. $\endgroup$ – jk_kim Sep 27 '17 at 23:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.