Commit 53b1a7a1 authored by Holger Brandl's avatar Holger Brandl

started pe-aligner wrapper

parent aec15c2a
#!/usr/bin/env kscript
//DEPS org.docopt:docopt:0.6.0-SNAPSHOT de.mpicbg.scicomp:joblist:0.6-SNAPSHOT de.mpicbg.scicomp.bioinfo:kutils:0.1-SNAPSHOT
// add docopts to local m2 index
// git clone && cd && mvn clean package install -Dmaven.test.skip=true
// kotlinc-jvm -classpath $(mvncp org.docopt:docopt:0.6.0-SNAPSHOT)
// update joblist
// cd ~/mnt/mack/bioinfo/holger/misc/Dropbox/cluster_sync/joblist/
import joblist.JobConfiguration
import joblist.JobList
import kutils.bashEval
import org.docopt.Docopt
import java.util.*
//val args = listOf("/projects/bioinfo/igenomes/Drosophila_melanogaster/Ensembl/BDGP6", "test.fastq.gz")
// basic usage tutorial
val usage = """
Use star to align fastq files against a genome
Usage: [options] <igenome> <fastq_files>...
--gtf <gtfFile> Custom gtf file instead of igenome bundled copy
--pc-only Use protein coding genes only for mapping and quantification
val doArgs = Docopt(usage).parse(args.toList()).map {
it.key.removePrefix("--").replace("[<>]", "") to it.value?.toString()
// extract all configuration parameters
val fastqFiles = (doArgs["fastq_files"] as String).split(" ").map { File(it) }
val igenome = File(doArgs["igenome"])
val star_index = File(igenome, "Sequence/StarIndex")
val gtfFile = if (doArgs["gtfFile"] != null) File(doArgs["gtfFile"]) else File(igenome, "Annotation/Genes/genes.gtf")
// check if gtf file exists
if (!gtfFile.isFile()) throw IllegalArgumentException("gtf '${gtfFile}' does not exist")
// make sure that STAR is in the PATH
if (bashEval("which STAR").exitCode != 0) throw IllegalArgumentException("gtf '${gtfFile}' does not exist")
// Check if STAR index is not present
if (!File("${star_index}/SA").isFile())
throw IllegalArgumentException("Missing STAR for ${star_index}/SA; use 'dge_create_star_index ${igenome}' to create it")
// dge_create_star_index ${igenome}'
println("running STAR using igenome '${igenome}' for the following files")
val jl = JobList(".starjobs")
//rm .starjobs
//fastqFiles=$(ls $baseDir/treps_pooled/*fastq.gz)
for (fastqFile in fastqFiles) {
println("submitting STAR job for $fastqFile")
// uniquely mapping reads only: http:/
//# tophat -p6 -G $gtfFile -g1 -o test $bowtie_gindex $fastqFile
// params:
// --outSAMstrandField intronMotif is required for cuffdiff compatiblity (xs flag)
// --outFilterIntronMotifs RemoveNoncanonicalUnannotated --outFilterType BySJout to get rid of artifical junctions
// --quantMode GeneCounts see!searchin/rna-star/GeneCounts/rna-star/gZRJx3ElRNo/p5FjBYKuY00J
// --outSJfilterCountUniqueMin see!topic/rna-star/_1BeAlGUmpA
val fastqBaseName =".fastq.gz", "")
val cmd = """
STAR --genomeDir $star_index --readFilesIn $fastqFile --runThreadN 6 --readFilesCommand zcat --outFileNamePrefix ${fastqBaseName}. --outSAMtype BAM SortedByCoordinate --outSAMstrandField intronMotif --sjdbGTFfile $gtfFile --outFilterIntronMotifs RemoveNoncanonicalUnannotated --outFilterType BySJout --quantMode GeneCounts --outFilterMultimapNmax 1 --outSJfilterCountUniqueMin 8 3 3 3
mv ${fastqBaseName}.Aligned.sortedByCoord.out.bam ${fastqBaseName}.bam
samtools index ${fastqBaseName}.bam
""".trimIndent(), "star__${fastqBaseName}", null, "medium", 5, "", better.files.File(File(".").toPath())))
//val requiresRerun: List<Job>? = jl.requiresRerun().
//if (requiresRerun.istEmpty()) throw RuntimeException("jobs failed: \n" +jl.printStatus())
//// estimate expresssion with
////Summarize multiple datasets at the same time:
//#featureCounts - t exon -g gene_id -a annotation.gtf -o counts.txt library1.bam library2.bam library3.bam
//// FeatureCounts --> Commented out because now part of STAR
//#bamFile = control_3.bam
//#featureCounts - t exon -g gene_id -a ${ gtfFile } - o feature_counts.txt ${ bamFile } - T 5
//#mysub "${project}__feature_counts" "featureCounts -t exon -g gene_id -a ${gtfFile} -o gene_counts.txt $(ls *bam) -T 5"-q medium | blockScript
// cleanup
bashEval("""rm - rf * STARgenome * .Log.progress.out _STARtmp * *Log.out""")
// Test for bam correlation (but don't wait for the results)
rend.R - e ${ System.getenv("NGS_TOOLS") } / dge_workflow / star_qc.R.
dge_bam_correlate &
// Condense counts into matrix
<project xmlns="" xmlns:xsi=""
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment