Implementera kvarvarande MEDUSA steg i StaG

[lis4matilda]

Implementera bowtie mapping, counting, combining, annotating

Medusa scripts:

• streamAligner.py
• combineCounts.py
• annotateCounts.py

input needed:

• filtered.humanremove.reads
• database index
• database annotation tsv file

Plan

• Starta med bowtie (Vi kanske kan implementera fler mappers senare).
  • Parat eller inte parat, det är frågan...
• Ge database index som input
• Ett folder-system per database (Ge folder samt filer namn som sätts i config-filen)
• Starta med den counting som finns i Medusa idag (lägg andra count-möjligheter som options)
• Statistics (maybe we can start with what is used in MEDUSA but if you want to adjust, please adjust.)
• Combine samples into one file (features som rader, prov som kolumner). Jag har skrivit hjälpskript till detta steget om ni vill kika.
• Summera statistik baserat på prov.
• Annotate. Collaps feature into the same annotation given in the annotation file.
  • input: columns to use for annotation (eg. 1,4,5).
  • Output: 1 combined countfile per annotation.

Filsystemet kommer bli (typ):

MEDUSA_databasnamn (folder)
       Counts (folder)
             counts
             counts.log
       Combine (folder)
            combine.feature
            combine.anno1
            combine.anno4
            combine.anno5
       summary.stats.all.samples

[boulund]

@lis4matilda, jag formaterade ditt inlägg lite så det blir lättare att läsa. Github stöder Markdown i inlägg, så det är enkelt att göra simpel formatering av text för att göra det tydligt att läsa. Hoppas jag fick formateringen rätt ungefär som du hade tänkt :).

Detta är en lite detalj, men jag funderar på om vi inte ska ha en output struktur för MEDUSA-grejerna som ser ut typ så här:

output_dir/
    MEDUSA/
        <db_name>/
            Counts/
                <sample>.counts.tab
                <sample>.counts.log
            Combine/
                <sample>.combine.feature
                <sample>.combine.anno1
                <sample>.combine.anno2
                ... osv...
        summary.stats.all.samples

Det blir nog lite enklare att ha en undermapp för varje databas man kört till MEDUSA, så ligger alla MEDUSA-output samlade i en mapp i den samlade outputmappen för hela workflowet.

Jag utgår ifrån att de MEDUSA-script som du nämner inte ändrats nämnvärt sen publikationen? Bifogar dem till det här inlägget för snabb access (men de hamnar nog förr eller senare i själva workflowet i någon form):

combineCounts.py.txt streamAligner.py.txt annotateCounts.py.txt

Gissningsvis kommer vi implementera bowtie2-mappningssteget utan streamAligner.py, för att hela det ramverket som finns där helt enkelt inte behövs när vi har Snakemake som ram runt hela vårt workflow. Sen kommer vi hugga relevanta delar ur combineCounts.py or annotateCounts.py.

[lis4matilda]

Det är helt ok :) Vi behöver inte kalla det Medusa alls. Kan bara inte komma på vad vi skall kalla det. Tänker att vi inte behöver implementera scripten alls, tror snakemake är bättre, bara att ha som referens om det behövs.

Lisa

5 apr. 2018 kl. 15:57 skrev Fredrik Boulund notifications@github.com:

@lis4matilda, jag formaterade ditt inlägg lite så det blir lättare att se. Github stöder Markdown annotering av poster, så det är enkelt att göra simpel formatering av text för att göra det tydligt att läsa. Hoppas jag fick formateringen rätt ungefär som du hade tänkt :).

Detta är en lite detalj, men jag funderar på om vi inte ska ha en output struktur för MEDUSA-grejerna som ser ut typ så här:

output_dir/ MEDUSA/

/ Counts/ .counts.tab .counts.log Combine/ combine.feature combine.anno1 combine.anno2 ... osv... summary.stats.all.samples Det blir nog lite enklare att ha en undermapp för varje databas man kört till MEDUSA, så ligger alla MEDUSA-output samlade i en mapp i den samlade outputmappen för hela workflowet. Jag utgår ifrån att de MEDUSA-script som du nämner inte ändrats nämnvärt sen publikationen? Bifogar dem till det här inlägget för snabb access (men de hamnar nog förr eller senare i själva workflowet i någon form): combineCounts.py.txt streamAligner.py.txt annotateCounts.py.txt Gissningsvis kommer vi implementera bowtie2-mappningssteget utan streamAligner.py, för att hela det ramverket som finns där helt enkelt inte behövs när vi har Snakemake som ram runt hela vårt workflow. Sen kommer vi hugga relevanta delar ur combineCounts.py or annotateCounts.py. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[boulund]

Exakt vad för output är det du förväntar dig från ett sånt här mappningssteg? Har för mig MEDUSA ger counts per annoterad region, stämmer det? I så fall kan vi ju bara köra pileup.sh från BBMap på outputfilen från bowtie2, så får du count/RPKM/osv per region i en fin textfil som är lätt att parsa. De flesta mappers klarar ju av att skapa SAM/BAM output, och om vi standardiserar output från alla mappningssteg som en pileup-output blir det ju enhetligt och fint och lättare att jobba downstream oavsett vilken mapper man använt för att mappa reads mot sin referensdatabas.

Det finns nu en branch i repot, bowtie2, som lagt in ett steg som mappar med bowtie2 till en databas som man specificerar i config.yaml. Om man anger en felaktig databas i configfilen startar inte pipelinen, utan den ger ett felmeddelande att den inte kunde hitta databasen.

[lis4matilda]

Hm, inte helt hundra på hur du menar, pileup när någon form av coverage om jag inte minns fel. Det viktigaste vi måste bestämma är hur vi assignar reads (best match, unique mapps, equal distribution). Tony tyckte vi kunde implementera alla, tror best match är implementerat i medusa nu. Sedan antar jag att det är som att mäta coverage per feature.

7 apr. 2018 kl. 13:18 skrev Fredrik Boulund notifications@github.com:

Exakt vad för output är det du förväntar dig från ett sånt här mappningssteg? Har för mig MEDUSA ger counts per annoterad region, stämmer det? I så fall kan vi ju bara köra pileup.sh från BBMap på outputfilen från bowtie2, så får du count/RPKM/osv per region i en fin textfil som är lätt att parsa. De flesta mappers klarar ju av att skapa SAM/BAM output, och om vi standardiserar output från alla mappningssteg som en pileup-output blir det ju enhetligt och fint och lättare att jobba downstream oavsett vilken mapper man använt för att mappa reads mot sin referensdatabas.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[boulund]

Jepp, du minns rätt: pileup från samtools beräknar coverage. pileup.sh från BBMap räknar ut mer än bara coverage. Den räknar dock ut statistik per referenssekvens, så om ni har komplicerade referensdatabaser med flera annoterade regioner per referenssekvens är det nog inte riktigt rätt verktyg. Det är sjukt smidigt i alla fall, och jag rekommenderar det starkt över samtools om du inte har nån annan preferens (eller nåt annat downstream-steg som absolut behöver output från samtools).

Det går redan nu att justera best match/unique mappings osv för bowtie2-mappningen med bowtie2:s egna inbyggda settings för det (bara att lägga till alla extra bowtie2-settings man vill ha med i en setting som heter "extra" i config.yaml). Notera dock att jag inte merge:at bowtie2-grenen med master-grejen ännu, så denna funktionalitet finns än så länge bara i bowtie2-grenen.

[lis4matilda]

Om vi har flera annoterade regioner per referenssekvens vet jag inte. Tror inte databaserna är så komplicerad. Fasta-filer med en annotering per sekvenssnutt. Skulle tro att det du tänker dig funkar bra. Du borde ha IGC-databasen, eller?

Vet inte om vi behöver alignments från mappningen, tror inte medusa sparar dem samt att de tar plats. Vad jag skulle gissa att man skulle vilja göra är att assembla de reads som mappar till specifika features men tror det är bättre att mappa om enbart på de features isf för att ta ut de reads än att alltid spara alla alignments.

7 apr. 2018 kl. 13:36 skrev Fredrik Boulund notifications@github.com:

Jepp, du minns rätt: pileup från samtools beräknar coverage. pileup.sh från BBMap räknar ut mer än bara coverage. Den räknar coverage per referenssekvens, så om ni har komplicerade referensdatabaser med flera annoterade regioner per referenssekvens är det nog inte riktigt rätt verktyg. Det är sjukt smidigt i alla fall, och jag rekommenderar det starkt över samtools om du inte har nån annan preferens (eller nåt annat downstream-steg som absolut behöver output från samtools).

Det går redan nu att justera best match/unique mappings osv för bowtie2-mappningen med bowtie2:s egna inbyggda settings för det (bara att lägga till alla extra bowtie2-settings man vill ha med i en setting som heter "extra" i config.yaml)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[boulund]

Det är helt upp till dig hur du vill göra med alla outputfiler från stegen. Givetvis kommer det ta en massa plats med en massa intermediära filer. Snakemake har inbyggda funktioner för att automatiskt ta bort intermediära filer när de inte längre behövs om man markerar dem som temp. Då håller Snakemake koll på om en fil fortfarande behövs av en downstreamprocess i workflowet, t.ex. en fastq med trimmade reads innan human-filtrering som man kanske inte behöver spara för all framtid (de går ju snabbt att göra om från rådatan om man verkligen behöver den), men den behövs ju fram tills man filtrerat bort human-DNA. När workflowet som helhet börjar bli mer komplett kan vi lätt se över vilka filer som faktiskt behöver sparas och vilka som kan raderas och justera det allt eftersom.

Jag tror du har helt rätt i att det är mycket effektivare att mappa mot en mindre databas med specifika features om man vill ha bara de reads som mappar till ett subset av feature-sekvenser. De flesta mappers kan väl mappa mot en referensdatabas och automatiskt producera en ny FASTQ-fil med bara de reads som mappade enligt matchningskriteria (BBMap gör det iaf).

[boulund]

Jag funderar på att splitta upp denna issue i ett par mindre bitar under "MEDUSA in Stag-mwc" milestonen. Denna kan få ligga kvar som en generell issue för allmän diskussion tills vidare.