text:
Nextflow (DSL v2): how to best synchronize multiple outputs from a single process
I have a workflow that needs to:

Generate .fastq files from .bams (while preserving the @RG group from the original bam)
Split the .fastqs
Align
etc.

The only way I could think to preserve the @RG group was to print it to a file during step 1 and emit the file. I can't figure out how to synchronize the emitted @RG file with the .fastqs, however. I've tried several variations without luck.
Questions:

How do I synchronize the emitted @RG file (sample_RG_file) with the split .fastqs?
Is there a better way?
Bonus question (may be worth a separate post): does Nextflow not stage files until the actual .bash script starts? I tried to operate on the sample_RG_file within the script using Groovy, but got an error that the file didn't exist. (The code was working before I found I needed to synchronize sample_RG_file with the

Other variations:

The code was 'working' before I found I needed to synchronize sample_RG_file with the split .fastqs. I was originally passing GENERATE_FASTQS_PROC.out.sample_RG as a separate parameter to ALIGN_FASTQ_BWA_PROC, but it was always passing the same (first) file.
I also tried GENERATE_FASTQS_PROC.out.sample_RG.flatten() in the .map operator, but that had the same error.

Current error:
The code below fails with the following error:

Error executing process > 'REALIGN_SAMPLES_WF:ALIGN_FASTQ_BWA_PROC (1)'

Caused by:
Not a valid path value type: groovyx.gpars.dataflow.DataflowBroadcast (DataflowBroadcast around DataflowStream[?])

Semi-pseudo code
Workflow
workflow my_WF {
    take:
        sample_input_file_ch
        original_ref

    main:
        GENERATE_FASTQS_PROC(sample_input_file_ch, original_ref)

        reads_per_job = 1_000_000
        fq_split_reads_and_RG_ch = GENERATE_FASTQS_PROC.out.fastqs
                        .splitFastq(by: reads_per_job, pe: true, file: true)
                        .map {
                            fq1, fq2 ->
                                tuple( fq1, fq2, GENERATE_FASTQS_PROC.out.sample_RG )
                        }
                        .view()

        ALIGN_FASTQ_BWA_PROC(fq_split_reads_and_RG_ch, ...)
}

GENERATE_FASTQS_PROC
process GENERATE_FASTQS_PROC {

    input:
        path(sample_input_file)
        val(original_ref)

    output:
        tuple path("**/*_R1.fastq"), path("**/*_R2.fastq"), emit: fastqs
        path("**/*_RG.txt"), emit: sample_RG

    script:
    """
    generate_fastqs.sh $sample_input_file $original_ref $task.cpus $mem_per_thread
    """
}

ALIGN_FASTQ_BWA_PROC
process ALIGN_FASTQ_BWA_PROC {

    /*
     * Publish results. 'mode: copy' will copy the files into the publishDir
     * rather than only making links.
     */
    // publishDir("${params.results_dir}/01-REALIGN_SAMPLES", mode: 'copy')

    label 'REALIGN_SAMPLES'

    input:
        tuple path(fq1), path(fq2), path(sample_RG_file)
        val(align_to_ref)
        val(align_to_ref_tag)
        val(original_ref)
        val(output_format)

    output:
        path '*.{bam,cram}', emit: final_alignment

    script:

    /*
     * Tried to read the RG file using Groovy, but
     * got error that the file didn't exist.
     *
     * Does Nextflow not stage files until the actual .bash script starts?
     */
//    println "sample RG file: $sample_RG_file"
//    file_lines = sample_RG_file.readLines()
//    sample_name = file_lines[0]
//    sample_RG = file_lines[1]

    """
    lines=()
    while read line
    do
        lines+=("\\\$line")
    done < $sample_RG_file

    sample_name="\\\${lines[0]}"
    sample_RG="\\\${lines[1]}"

    bash realign_bwa-split.sh $fq1 $fq2 \\\$sample_name \"\\\$sample_RG\" ${align_to_ref} ${output_format} $task.cpus
    """
}

Update
Per @Steve's suggestion, I'm trying to have GENERATE_FASTQS_PROC create the split .fastq files using Unix split, but I'm still struggling with various concepts.
Supposed I have this emitted from my process:
[sample_A_RG.txt, [s1_R1.split_00.fq, s1_R2.split_00.fq, ..., s1_R1.split_XX.fq, s1_R2.split_XX.fq]]`

In my mind, I want to pull out s1_RG.txt file from the Channel and extract the @RG. Then I want to call:
[<all_fastqs_for_sample_A>]
    .map{
        // pull out split ID numbers
    }
    .groupTuple()

I believe this would group _R1 and _R2 for each split set, and then call
[[00.fq1, 00.fq2], ..., [XX.fq1, XX.fq2]]
    .combine(sample_RG_var)

to get [[@RG, 00.fq1, 00.fq2], ..., [@RG, XX.fq1, XX.fq2]].
Seems like it would be so much easier if I could use Channel factories and operators (e.g., .fromFilePairs().combine( rg )) in the output:.
Getting files into a process is relatively straightforward, but getting them out in a reasonable manner is less so.

Nextflow channels guarantee that items are delivered in the same order as they are sent. So if a process declares two (or more) output channels, the items that are emitted should already be synchronized. This is true unless of course one of the channels sets the optional output attribute.

How do I synchronize the emitted @RG file (sample_RG_file) with the
split .fastqs?

Since the splitFastq operator can create n outputs for each input, you need some way to combine the @RG files using some key. This key, for example, could be the basename of the sample input file. This, I think, would work nicely given the input file should just have a '.bam' extension. Another option, would be to supply some other value in the input declaration and just use that in your output declaration. If we went with the former, your output declaration might look like:
output:
tuple val(sample_input_file.baseName), path("sample_R1.fastq"), path("sample_R2.fastq"), \
    emit: fastqs
tuple val(sample_input_file.baseName), path("sample_RG.txt"), \
    emit: sample_RG

Then, the body of your workflow might look like:
workflow {

    reads = GENERATE_FASTQS_PROC(sample_input_file_ch, original_ref)

    reads.fastqs
        .splitFastq( by: params.reads_per_job, pe: true, file: true )
        .combine( reads.sample_RG, by: 0 )
        .map { sample, fq1, fq2, rg_txt -> tuple(fq1, fq2, rg_txt) }
        .set { fq_split_reads_and_RG_ch }

    ...
}

You'll notice that I removed your fuzzy glob patterns from your output declaration above. Your workflow definition suggests that the generate_fastqs.sh script produces a single pair of FASTQ files and a single readgroup file for a given set of inputs. However, as you know, glob patterns can match one or more output files. When declaring output files, it's always best to name these guys explicitly. If you do use a glob pattern and for some reason your process was to produce two files that match a given pattern, you'll instead get back a tuple, which would break the pipeline.

Is there a better way?

Since you need the @RG tags with the split FASTQ files, I don't think there is. However, I generally avoid using the splitting operators on large files. My preference is always to pass "work" off (including IO) to a separate process. A separate process can also provide increased flexibility if/when requirements change.

Does Nextflow not stage files until the actual .bash script starts? I tried to operate
on the sample_RG_file within the script using Groovy, but got an error
that the file didn't exist.

Nextflow evaluates the code in the script definition before returning the script string to be executed. This is why it's possible to implement a conditional script. Basically, you want to avoid all IO operations here. The script string to be executed has not even been decided at this point.

Update
I'm not sure how you plan to operate per sample. IIUC, [<all_fastqs_for_sample_A>] is basically a Channel of lists, where each list contains all split FASTQ files for a given sample. Using map here to pull out the split ID numbers followed by the groupTuple operator might not give the expected results. You might also run into the problems with calling groupTuple, since, without a size attribute, it will basically wait for all outputs before emitting the grouped result:

You should always specify the number of expected elements in each
tuple using the size attribute to allow the groupTuple operator to
stream the collected values as soon as possible. However there are use
cases in which each tuple has a different size depending on the
grouping key. In this case use the built-in function groupKey that
allows you to create a special grouping key object to which it’s
possible to associate the group size for a given key.

Part of the problem is that we don't know what the generate_fastqs.sh script does and how the other scripts work. If you are aligning using BWA MEM and have samtools available in your $PATH, it might even be easier to produce interleaved output (since BWA MEM can read interleaved FASTQs) using the samtools fastq command. This would save some work having to group each list of FASTQs into pairs, ensuring they are ordered as expected. This would be reasonably involved.
Not sure what BAMs you need to process, or if you're looking for a more generic solution, but I thought I'd throw in an alternative solution you might find works well for you or that you may not have yet considered. Since BAMs often contain multiple readgroups, I find that if I split by readgroup (using samtools split, flatten the output Channel, then samtools collate followed by samtools fastq, I can achieve somewhat uniform mapping walltimes (assuming the input BAMs are similar) without needing to chunk the reads and then merge the alignments back together. If your input BAMs are not similar at all, then you will likely need to chunk to generate uniform alignment jobs.
Below is an example, that shows how to create a groupKey, followed by the transpose, join/combine and groupTuple pattern. There are some caveats and there is some room for optimization:
nextflow.enable.dsl=2

params.original_ref = './ref.fasta'
params.reads_per_job = 10000

process samtools_view_header {

    tag { "${sample_name}:${bam.name}" }

    cpus 1
    memory 1.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.header.txt")

    """
    samtools view \\
        -H \\
        -o "${bam.baseName}.header.txt" \\
        "${bam}"
    """
}

process samtools_name_sort {

    tag { "${sample_name}:${bam.name}" }

    cpus 4
    memory 8.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.nsorted.bam")

    script:
    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    samtools sort \\
        -@ "${task.cpus - 1}" \\
        ${mem_per_thread} \\
        -n \\
        -o "${bam.baseName}.nsorted.bam" \\
        -T "${bam.baseName}.nsorted" \\
        "${bam}"
    """
}

process samtools_fastq {

    tag { "${sample_name}:${bam.name}" }

    cpus 1
    memory 2.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.fastq")

    """
    samtools fastq \\
        -O \\
        -T RG,BC \\
        -0 /dev/null \\
        "${bam}" \\
        > \\
        "${bam.baseName}.fastq"
    """
}

process split_fastq {

    tag { "${sample_name}:${fastq.name}" }

    cpus 1
    memory 1.GB

    input:
    tuple val(sample_name), path(fastq)

    output:
    tuple val(sample_name), path("${fastq.baseName}.${/[0-9]/*5}.fastq")

    """
    split \\
        --suffix-length=5 \\
        --additional-suffix=".fastq" \\
        -d \\
        -l "${params.reads_per_job*4}" \\
        "${fastq}" \\
        "${fastq.baseName}."
    """
}

process bwa_index {

    tag { fasta.name }

    cpus 1
    memory 12.GB

    input:
    path fasta

    output:
    tuple val(fasta.name), path("${fasta}.{amb,ann,bwt,pac,sa}")

    """
    bwa index "${fasta}"
    """
}

process bwa_mem {

    tag { "${sample_name}:${fastq.name}" }

    cpus 8
    memory 12.GB

    input:
    tuple val(sample_name), path(fastq), path(header)
    tuple val(idxbase), path(bwa_index)

    output:
    tuple val(sample_name), path("${fastq.baseName}.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    bwa mem \\
        -p \\
        -t ${task_cpus} \\
        -C \\
        -H <(grep "^@RG" "${header}") \\
        "${idxbase}" \\
        "${fastq}" |
    samtools view \\
        -1 \\
        -o "${fastq.baseName}.bam" \\
        -
    """
}

process samtools_coordinate_sort {

    tag { "${sample_name}:${bam.name}" }

    cpus 4
    memory 8.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.csorted.bam")

    script:
    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    samtools sort \\
        -@ "${task.cpus - 1}" \\
        ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
        -T "${bam.baseName}.csorted" \\
        "${bam}"
    """
}

process samtools_merge {

    tag { sample_name }

    cpus 1
    memory 4.GB

    input:
    tuple val(sample_name), path(bams)

    output:
    tuple val(sample_name), path("${sample_name}.bam{,.bai}")

    """
    samtools merge \\
        -c \\
        -p \\
        "${sample_name}.bam" \\
        ${bams}
    samtools index \\
        "${sample_name}.bam"
    """
}

workflow {

    original_ref = file( params.original_ref )
    bwa_index( original_ref )

    Channel.fromPath('./data/*.bam')
        | map { tuple( it.baseName, it ) }
        | set { sample_bam_files }

    samtools_view_header( sample_bam_files )

    sample_bam_files
        | samtools_name_sort
        | samtools_fastq
        | split_fastq
        | map { sample_name, fastq_files ->
            tuple( groupKey(sample_name, fastq_files.size()), fastq_files )
        }
        | transpose()
        | combine( samtools_view_header.out, by: 0 )
        | set { realignment_inputs }

    bwa_mem( realignment_inputs, bwa_index.out )
        | samtools_coordinate_sort
        | groupTuple()
        | map { sample_name, bam_files ->
            tuple( sample_name.toString(), bam_files )
        }
        | samtools_merge
        | view()
}

There are some caveats however, in no particular order:

It unravels your shell scripts. I think this is a better solution though because it gives you control over each of the components.

No CRAM support. But this should be easy to add if needed.

A single process to read a BAM header is wasteful, especially if the input BAM is large. More time will be spent localizing the file than actually doing work.

Replace samtools sort -n with samtools collate. There's really no need to name sort. You may prefer to pipe the output of samtools collate into samtools fastq for avoid writing an extra file. This might be useful if you need to operate on an extremely large number of large BAMs.

The 'samtools_fastq' step should write gzipped output with low compression. This means replacing coreutils split with something else. A solution using Biopython could be useful here to read and write gzipped output in a single pass.

The last 'map' command in the workflow is unnecessary. I just thought I'd show that calling .toString() on the GroupKey gives you back an actual String object.

Also, it's not necessary to pass around a tuple of 'sample name' and paths as process inputs and outputs. There's no harm in doing this though. When you are more confident with Nextflow Channels, refactoring these processes should be easy.

My samtools sort processes ensure that they do not oversubscribe the nodes that they land on. They use approximately 320% - 380% CPU, depending on the size of the BAM. The resource requests for the other processes are not optimized at all.

With newer versions of samtools, BAM indexing can be done at the same time as merging, with the --write-index option.

Thanks, @Steve. You've made several important points. Question about your suggested workflow code: wouldn't .splitFastq fail since reads.fastqs contains a tuple including the sample_file_input.baseName? I've been battling how to work with the tuples of tuples in various scenarios (e.g., in this situation and after a .groupTuple()).

Also, per your suggestion, I've been working to have GENERATE_FASTQS_PROC split the .fastqs itself and emit those. That's causing me a headache. The level of abstraction in Nextflow has really made these details quite elusive for me.

I added an update to my question explaining what I feel like I should be able to do and what I ultimately want to get.

Thanks, @Steve. I think the interleaved option might be the simplest. Great suggestion. Still, I'm shocked by how difficult it is to get files out of a process without mixing them between samples. What kills me is, I know all files for a single sample are safely contained in a unique directory, but I'm forced to mix them up with files from other samples in the same Channel, so I'm left with convoluted methods to try to distinguish them from each other. If I emit the path to the directories and then use Channel.fromPath(), it would be sooooo easy. But that's not possible, AFAIK.

generate_fastqs.sh takes a single .bam -> sorts by name -> converts to .fastq (currently w/ bedtools), and splits them with Unix split. The files are stored safe and sound, all tidy-like in a unique directory, where I have clear control over what's generated and mixed together (unlike Nextflow Channels). Forgive the frustration. :-)

Sorry @MarkEbbert, I was meant to let you know yesterday that I had updated my answer but got pulled in another direction before I could finish the update. WRT your first question above, splitFastq will work fine with a tuple that includes a sample identifier. For example, it will work with the output of fromFilePairs with flat: true.

@MarkEbbert I struggled with the same thing when I first started with NF. I found that outputting tuples of files, prefixed with some identifier (usually a sample name) helps considerably when working with multiple channels.

This was a huge help, @Steve. You've also made several additional helpful suggestions. It had never dawned on me that samtools was oversubscribing because -@ is for additional threads, though that should have been completely obvious. Really appreciate your help. I'm still not a fan of NF output Channels, but maybe I'll get there. :-)

Update 3
While the code in Update 2 works, @Steve pointed out that my approach may not be ideal in cloud environments (see comments). He also provided a Nextflow solution that works. I'm still not a fan of how Nextflow output: Channels work, but that may be from my inexperience.
Update 2
Pretty sure I got it this time. The following code works. I'm getting all of the tuples I was expecting.
Here's the updated code for my workflow to create the tuples from each line in the sample_tuple_file I created from a bash for loop. I spent quite a while trying to get a combination of .splitText() and .splitCSV() to work, but for some reason, you cannot chain as follows: it.splitText().splitCSV().
     GENERATE_FASTQS_PROC.out.sample_tuple_file
         .map {
            it.splitText()
         }
         .flatten()
         .map {
             toks = it.trim().split(',')
             sample_name = toks[0]
             fq = toks[1]
             rg = toks[2]
             tuple(sample_name, fq, rg)
             }
         .view()
        .set { sample_fq_tuples_ch }

Update 1
I thought I had it, but I'm actually only getting ONE of the tuples from my GENERATE_FASTQS_PROC.out.sample_tuple_file.map { } statement (the last one). And I have no idea why...

Original posted answer
Despite @Steve's exceptional help and patience with me, I have really battled to figure out how to work with the Nextflow output: Channels, but simply could not figure out any solution using them.
I finally figured out a super simple method that is much easier than dealing with Nextflow Channels, which I share here in case it will help others. I personally will probably use this method from now on, unless someone can convince me that's a bad idea. In my view, using my approach gives me complete control over the way outputs and inputs are handled, and I don't have to worry whether Nextflow is doing what I hope it's doing.
generate_fastqs.sh
generate_fastqs.sh takes a single .bam -> sorts by name -> converts to .fastq (now w/ samtools fastq, thanks to @Steve), and splits them with Unix split. Then, at the end of the script, I create a single file with all of the tuples that I want to pass to my next process.
Here are the key lines for these steps (including exit codes to ensure Nextflow receives error notifications). The last step is the key (for me, anyway):
...

#####################
# Sort .bam by name #
#####################
if ! time samtools sort -@ $n_threads -n -m ${mem_per_thread}G $sample_input > $sample_sorted_by_name; then
    echo "ERROR: sort failed"
    exit 1
fi

...

#############################
# Create interleaved .fastq #
#############################
if ! time samtools fastq -o $fq -s /dev/null -0 /dev/null --reference $original_ref $sample_sorted_by_name; then
    echo "ERROR: bamtofastq failed"
    exit 1
fi

...

################
# Split fastqs #
################
if ! time split --additional-suffix=".fastq" -dl $n_lines $fq "${split_dir}/${sample_name}.interleaved_R1_R2.split."; then
    echo "ERROR: split failed"
    exit 1
fi

...

###########################################################
# Store sample name, .fastq, and RG to file for alignment #
###########################################################
RG_file="${sample_name}.tuples.txt"

# Create empty file
> $RG_file
for fq_file in $PWD/${split_dir}/*
do
echo "${sample_name},$fq_file,$RG" >> $RG_file    # <-- Creating future tuples in file, here (comma separated).
done

GENERATE_FASTQS_PROC
In GENERATE_FASTQS_PROC, emit one simple file--the tuple file I created myself with a simple bash loop (so easy and feels so good):
    emit:
        path("*.tuples.txt"), emit: sample_tuple_file

Workflow
In my workflow, I call GENERATE_FASTQS_PROC, access the file from the Channel, read each line, and create a Nextflow tuple for each line (again, the lines I so easily created in a bash for loop).
    /*
     * Generate fastq files from the original .(cr|b)am file and split into
     * mini .fastq files with `n_reads_per_run` in each mini .fastq
     */
    GENERATE_FASTQS_PROC(sample_input_file_ch, n_reads_per_run, original_ref)

    /*
     * Get sample name, fastq, and read group (RG) from the tuple file using Nextflow
     */

     GENERATE_FASTQS_PROC.out.sample_tuple_file
        .map{
            sample_tuple_file ->
                sample_tuple_file.withReader {          // Open file
                    it.eachLine {                       // Loop over each line in file
                        line ->
                            toks = line.split(',')      // Split by comma
                            sample_name = toks[0]       // Get each item
                            fq = toks[1]
                            rg = toks[2]

                            tuple(sample_name, fq, rg)  // Create the simple tuple that has eluded me for two days!
                    }
                }
        }
        .set { sample_fq_tuples_ch }                    // Assign it to a new Channel name and move along!

    /*
     * Align fastq file pairs
     */

    ALIGN_FASTQ_BWA_PROC(sample_fq_tuples_ch, align_to_ref,
                             align_to_ref_tag, original_ref, output_format)
```

Ugh. I thought I had it, but I confused myself with my own print statements. I'm actually only getting ONE of the tuples from my GENERATE_FASTQS_PROC.out.sample_tuple_file.map { } statement. And I have no idea why...

Pretty sure I got it working with update #2. Hopefully others can learn from my mistakes and final victory.

I can see you have put in quite a bit of effort here, but this is not something I would recommend to others. This is a solution that will work on a local filesystem only, where all files exist under a common root. This is what makes it possible to pick them up later on if you have a pointer to them (your tuples.txt contains file pointers). Since the required files are not declared in the output block, Nextflow will not attempt to "unstage" them (upload them to a remote object store, like a bucket) when you deploy your workflow to the cloud. They will just be lost.

Thanks for including your generate_fastqs.sh. I have a better picture now of what you are trying to do. I will try to update my answer with a working solution soon. I just need to grab a couple of test BAMs. I should have a couple somewhere..

True. I hadn’t considered cloud file systems. I was worried you’d burst my bubble. :-)

Couldn’t I still emit the files you mention, but use my method to generate tuples, etc? Or, couldn’t I publish results and access them that way? I realize this breaks convention, but it appears that, while Nextflow has many advantages,  Nextflow’s output Channels are fundamentally convoluted, flawed, and maybe even dangerous. I’d love to be proven wrong. I don’t know Nextflow’s internals, but seems like this could be easily fixed. e.g., if every process instance assigned an ID to the process that we could filter by upon return, I could trust that NF is giving me what I think I’m getting.

You could perhaps. But that would not be a pleasant experience at all and would produce a convoluted solution. I don't have any real problem with Nextflow Channels. I think they're great and they make writing complex workflows possible. I've updated my answer above. I think it's ok as an example but there are some caveats listed. You may or may not need to address these. But it will likely need some tinkering before you can run with it.


text:
Nextflow: dynamically setting resource directives for a process
I'm trying to test variations of CPU and memory for a specific process (samtools sort) to measure computational efficiency (measured only by time). To do this, I'd like to set the CPU and memory for a Nextflow process dynamically.
Ultimate question: How can I dynamically set CPU and memory for many instances of a single Nextflow process definition?
Attempt 1
The most obvious solution was to submit distinct Nextflow workflows for each combination of CPU and memory by specifying the CPU and memory on the command line (i.e., when calling nextflow run ...), where my Nextflow workflow has global parameters params.samtools_total_threads and params.samtools_mem_per_thread.
This didn't work perfectly because there appeared to be 'collisions' between parent Nextflow instances trying to lock the history file. I got the following error:
Can't lock file: /path/to/project/.nextflow/history -- Nextflow needs to run in a file system that supports file locks

I could try delaying submissions, but I suspect that will still result in collisions.
Here's the SLURM script that I'm using:
TOTAL_THREADS=$1
MEM_PER_THREAD=$2

nextflow run 02-SAMTOOLS_SORT.nf \
    -with-report "reports/report.${TOTAL_THREADS}_threads.${MEM_PER_THREAD}_GB_per_thread.html" \
    -with-trace "traces/trace.${TOTAL_THREADS}_threads.${MEM_PER_THREAD}_GB_per_thread.txt" \
    --samtools_total_threads $TOTAL_THREADS \
    --samtools_mem_per_thread $MEM_PER_THREAD

My workflow had the following global parameters:
/*
 * Make this pipeline a nextflow 2 implementation
 */
nextflow.enable.dsl=2

...

params.samtools_total_threads = '4'
params.samtools_mem_per_thread = '8'

...

Attempt 2
After Attempt 1 failed, I figured the most 'Nextflow-like way' would be to generate tuples for each combination of CPU and memory that I want to test. Unfortunately, it doesn't appear that it's possible to set CPU and memory directives using parameters passed into the process, but only with global variables.
I get the following No such variable error with the below code (note: most of this code is courtesy of @Steve):
N E X T F L O W  ~  version 21.04.3
Launching `02-SAMTOOLS_SORT.nf` [desperate_aryabhata] - revision: fb9e79055a
No such variable: total_threads

 -- Check script '02-SAMTOOLS_SORT.nf' at line: 33 or see '.nextflow.log' file for more details

Process code
process samtools_coordinate_sort_proc {

    tag { "${bam.baseName}" }

    executor='slurm'
    queue = 'normal'
    cpus = total_threads
    memory = total_mem.GB
    clusterOptions = "--time 5:00:00 --account <account>"

    /* delete files upon completion (I think) */
    cleanup = true

    input:
    path(bam)
    val(total_threads)
    val(total_mem)

    output:
    path("${bam.baseName}.csorted.bam")

    script:

    def additional_threads = task.cpus - 1

    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    println "avail_mem: $avail_mem"
    def mem_per_thread = avail_mem ? "${avail_mem}G" : ''
    println "mem_per_thread: $mem_per_thread"

    """
    samtools sort \\
        -@ "${additional_threads}" \\
    -m ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
    -T "${bam.baseName}.csorted" \\
        --write-index \\
        "${bam}"
    """
}

Possible attempt 3
Of course, it would be super easy to do this with a couple of bash scripts, but I'm trying to do it in Nextflow for two reasons:

To learn Nextflow better
Because the traces provided by Nextflow are super convenient (though I could get this information from SLURM, too).

It's possible to specify dynamic computing resources using a closure. For example, this sets the cpus and memory directives using the values supplied via the input declaration:
nextflow.enable.dsl=2

process test_samtools {

    echo true

    cpus { cpus }
    memory { memory }

    input:
    tuple path(bam), val(cpus), val(memory)

    script:
    def additional_threads = task.cpus - 1

    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    echo samtools sort \\
        -@ "${additional_threads}" \\
        ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
        -T "${bam.baseName}.csorted" \\
        --write-index \\
        "${bam}"
    """
}

workflow {

    cpu_ch = Channel.of( 1, 2, 3, 4, 5)
    memory_ch = Channel.of( 1.GB, 2.GB, 4.GB )

    Channel.fromPath( './data/*.bam' )
        | combine( cpu_ch )
        | combine( memory_ch )
        | test_samtools
}

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Nextflow scheduler job doesn't end upon successful workflow completion
Bug fixed
This issue appears to have been fixed in latest Nextflow: https://github.com/nextflow-io/nextflow/issues/2693

I've been using Nextflow for various workflows on a cluster using the SLURM scheduler. I've noticed that the 'parent' job does not shut down when the Nextflow workflow finishes successfully. I can cancel the job, but then Nextflow doesn't get a chance to clean things up and also doesn't write out the trace, report, etc. The parent job does end properly if the workflow ends with an error.
Question: Why doesn't the parent Nextflow job end upon successful completion?
I saw a post on the Nextflow GitHub where the original poster describes various acrobatics to get the parent job to exit cleanly, but this is something that Nextflow should handle.
Here's my SLURM submission script that I submit using sbatch run_workflow.slurm
run_workflow.slurm
#!/bin/bash
#SBATCH --time=60:15:00                             # Time limit for the job (REQUIRED).
#SBATCH --job-name=NF_Parent_Identify_D&C_regions   # Job name
#SBATCH --ntasks=1                                  # Number of cores for the job. Same as SBATCH
#SBATCH --mem=5G                                    # Total memory requested
#SBATCH --partition=<partition>                     # Partition/queue to run the job in. (REQUIRED
#SBATCH -e slurm-%j.err                             # Error file for this job.
#SBATCH -o slurm-%j.out                             # Output file for this job.
#SBATCH -A <account>                                # Project allocation account name (REQUIRED)

module load ccs/java/jdk1.8.0_202
nextflow run identify_dark_and_camouflaged_regions.nf \
    -with-report \
    -with-trace \
    -with-dag idcr-workflow.png \
    -with-timeline \
    -resume

Update 1
I'm adding example Nextflow output from a different workflow (I don't have the SLURM output from the one I mention above) where the workflow has completed, but the parent job is still running. In this case, there were errors, but I added errorStrategy 'ignore' to all of the processes listed:
executor >  slurm (876)
[06/c5672b] process > samtools_csort_proc (A-CUHS... [100%] 304 of 304, faile...
[45/26a046] process > Lsamtools_csort_proc (A-CUH... [100%] 304 of 304, faile...
[b5/0f2d7c] process > sambamba_csort_proc (A-CUHS... [100%] 304 of 304, faile...

So, regardless of whether the workflow completed successfully, Nextflow has nothing left to do, but it's still running.
Update 2 (adding partial .nextflow.log file)
Per @Steve's point, I checked the .nextflow.log file and am including part of the file, here. In this case, I don't think Nextflow has any tasks related to publishDir because I'm not using that in this workflow. The job that's still running is from the Nextflow workflow in my GitHub repository for my answer posted in my related question Samtools sort: most efficient memory and thread settings for many samples on a cluster. Whether it's reproducible, I haven't checked. Not sure how common this issue is.
In either case, my interpretation of the .nextflow.log file below is that Nextflow thinks it's hung, too. I just don't know why, or how to exit cleanly.
Partial view of .nextflow.log
...

Feb-28 05:29:04.607 [Task monitor] DEBUG n.processor.TaskPollingMonitor - !! executor slurm > tasks to be completed: 1 -- submitted tasks are shown below
~> TaskHandler[jobId: 662274; id: 876; name: sambamba_csort_proc (A-CUHS-CU003128-BL-COL-49696BL1.unsorted;CPUs:2;MEM_PER_THREAD:100GB); status: RUNNING; exit: -; error: -; workDir: /path/to/Samtools_sort_optimization_test/SORT_COMPARISON/work/b5/0f2d7c5d7f81a6e62dedb08bf6694a started: 1646038174334; exited: -; ]
Feb-28 05:30:54.613 [Task monitor] DEBUG n.processor.TaskPollingMonitor - Task completed > TaskHandler[jobId: 662274; id: 876; name: sambamba_csort_proc (A-CUHS-CU003128-BL-COL-49696BL1.unsorted;CPUs:2;MEM_PER_THREAD:100GB); status: COMPLETED; exit: 0; error: -; workDir: /path/to/Samtools_sort_optimization_test/SORT_COMPARISON/work/b5/0f2d7c5d7f81a6e62dedb08bf6694a started: 1646038174334; exited: 2022-02-28T10:30:53.516015Z; ]
Feb-28 05:34:04.641 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:39:04.647 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:44:04.653 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:49:04.658 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:54:04.662 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled

< above line repeats hundreds of times >

...

Can you please create an answer with your accepted solution, so that the StackExchange algorithms don't confuse this with an unanswered question?

How do you know your workflow actually finishes successfully? One way to check the workflow termination status is to use the workflow.onComplete and workflow.onError event handlers. Consider resuming your workflow with something like the following appended to your script:
workflow.onComplete {

  log.info "Workflow completed at: ${workflow.complete}"
  log.info "Time taken: ${workflow.duration}"
  log.info "Execution status: ${workflow.success ? 'success' : 'failed'}"
}

workflow.onError {

  log.info "Execution halted: ${workflow.errorMessage}"
}

Note that the workflow event handlers may also be defined inside your nextflow.config if preferred.

Thanks, @Steve. I added an example to my post. Seems that one way or another, the example I show above should have completed because Nextflow has nothing left to do.

@MarkEbbert It's hard to say what the problem is here exactly without the nextflow.log or a reproducible example. But just because all of the tasks have completed (regardless of the exit status) doesn't necessarily mean Nextflow has nothing left to do. For example, the publishDir directive works asynchronously and can be used (as you know) to 'copy' task outputs. If these are large and there are many of these, Nextflow will hang around copying files from the queue until well after the last task has completed.

Thanks, @Steve. I just added part of the .nextflow.log file. Looks like Nextflow thinks it's hung, too? I haven't tested whether this is reproducible, by this is the code for my related question. I included a link to it above. This one may just remain a mystery.

I can't edit my above comment, but should read "..., but this is from the code..."

@MarkEbbert Correct - this looks like a bug. It might have something to do with the errorStrategy set to 'ignore' for all processes, potentially. Basically, you need to hit this condition to break out of the poll loop. From the partial log you've posted, it looks like the runningQueue is empty, but I don't know about the pendingQueue.size() or if session.isTerminated().

@MarkEbbert You might like to open an issue with the full log attached.

Thanks, @Steve. Submitted: https://github.com/nextflow-io/nextflow/issues/2693. I'll report back if I learn something.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
hg38 GTF file with RefSeq annotations
I'm not sure what I'm missing, but I'm struggling to find an official hg38 GTF file with RefSeq annotations. I'd like to provide the GTF to Salmon to get gene-level annotations.
Here's Salmon's help info for --geneMap:

File containing a mapping of
                                          transcripts to genes.  If this file is
                                          provided Salmon will output both
                                          quant.sf and quant.genes.sf files,
                                          where the latter contains aggregated
                                          gene-level abundance estimates.  The
                                          transcript to gene mapping should be
                                          provided as either a GTF file, or a in
                                          a simple tab-delimited format where
                                          each line contains the name of a
                                          transcript and the gene to which it
                                          belongs separated by a tab.  The
                                          extension of the file is used to
                                          determine how the file should be
                                          parsed.  Files ending in '.gtf', '.gff'
                                          or '.gff3' are assumed to be in GTF
                                          format; files with any other extension
                                          are assumed to be in the simple format.
                                          In GTF / GFF format, the
                                          "transcript_id" is assumed to contain
                                          the transcript identifier and the
                                          "gene_id" is assumed to contain the
                                          corresponding gene identifier.

Here's what I've tried:

NCBI seems to only provide GFF files. I know you can convert between the two, but that hasn't worked well for me in the past. Both file formats allow a lot of freedom, which makes conversions sloppy. Salmon, specifically, is looking for the gene_id and transcript_id annotations.
UCSC provides GTF files from RefSeq, but the gene_id annotation is identical to the transcript_id annotation (i.e., it's the NM number). Or maybe there's an option I'm missing.
Ensembl provides exactly the format and all of the information I could hope for, except the transcript_id is the Ensembl ID (ENST#), not RefSeq (NM#).

Can anyone tell me how to get a GTF file with hg38 RefSeq annotations?

refseq IDs are per-transcript, so it's expected that the gene and transcript IDs will be the same. Given that, what would you like as the gene_id?

@Devon Ryan I'm guessing the EntrezGene ID or HUGO symbol?

@Mark Ebbert: What are your down stream applications? The as far as I'm aware, the salmon output doesn't include the gene_id, so it shouldn't matter that it is the same as the transcript_id.

I appreciate your help. I added Salmon's help on --geneMap. They expect gene_id to contain a "gene identifier". I understand that could be several things, but I would prefer the HUGO symbol, as @IanSudbery suggested. In my perusing, I thing Salmon should probably target the gene_name, at least how Ensembl formats the GTF. In either case, I'm surprised by @DevonRyan's comment that the gene_id would equal the transcript_id. Not arguing that's the case, but what good does that do anyone?? :-)

@IanSudbery, I ultimately want to aggregate transcript expression into the parent genes. I found last night that one of the Salmon authors now recommend tximport (https://github.com/COMBINE-lab/salmon/issues/98), which I will use. But this isn't the first time I've "needed" a RefSeq GTF and couldn't find one. Just surprises me, but maybe that indicates I need to assess my approach?

You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format
The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our internally stored genePred format into GTF/GFF3.
This is the command I usually give when we get this question on our mailing list:

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout

Change stdout to the output filename you want in the last command to get an hg19 refGene GTF file:

chr1    hg19.refGene.ucsc   transcript  11869   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    11869   12227   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    12613   12721   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    13221   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   transcript  11874   14409   .   +   .   gene_id "DDX11L1"; transcript_id "NR_046018";  gene_name "DDX11L1";
...

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

General questions: genome@soe.ucsc.edu
Questions involving private data: genome-www@soe.ucsc.edu
Questions involving mirror sites: genome-mirror@ose.ucsc.edu

It would be rather nice if the "GTF" format output from the UCSC table browser simply returned exactly this (including the gene entries, which are always missing). This is the sort of thing people are actually expecting from it.

I agree that this would be ideal. The problem is that the gene_id's returned may not be what people want. For example, if you were trying to get a GTF output for the wgEncodeGencodeBasicV27 table, would you want a gene_id like DDX11L1 or a gene_id like ENSG*? The other issue is the plain genePred table format doesn't even have a name2 field, so what do we assign as gene_id in this case? What if the table is some other variation of genePred that also doesn't have a name2? We use many different table variations, and it is hard to support a GTF output for all of them.

@ChristopherLee, that's perfect, thank you. Regarding your response to @DevonRyan's comment, I can appreciate that it's difficult (and maybe unrealistic) to support all possible variations, but I'd also argue that HUGO is probably what most people want. In either case, I'm certain that reporting the transcript ID for both gene_id and transcript_id is not ideal! :-)

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Which quality score encoding does PacBio use?
Do you know which quality score encoding PacBio uses now? I know some of their file formats have changed in the past year or two, but I haven't found much on their quality score encoding.
The most recent answer I found is from 2012, where one user states it's "phred style", but doesn't state whether it's 33 or 64: DeNovo assembly using pacBio data

PacBio does use PHRED 33, but it turns out the question may be irrelevant for the newer PacBio Sequel Sequencer, because it reports all base qualities as PHRED 0 (ASCII !). The RS-II reports PHRED 33 quality scores.
The Sequel provides data in the form of subreads, which are the circular consensus sequences (CCS) from a single zero-mode waveguide (ZMW). Essentially, a single circularized DNA molecule enters a ZMW well, and is sequenced as many times as possible. (Sometimes multiple molecules get into a single well, but that's not the normal case.) The number of times this single molecule can be sequenced depends on several variables, including molecule length. The consensus sequence for this single molecule is then reported out as a subread.
As I started playing with the quality scores in my data, I realized the quality score for every base is !, which equates to a PHRED quality of 0. I contacted someone at PacBio who informed me that they learned the base quality values were not any more useful than the bases themselves for generating a high-quality circular consensus sequence. Thus, rather than spend the computational effort to calculate the scores, they skipped it, and reported ! for every base. They also report an overall read quality of 0.8, but this is simply a placeholder. They are working to make this information more readily available.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Estimate the length of poly-A tails from randomly-primed RNAseq data
So a poly-A tail is a long chain of adenine nucleotides that is added to a messenger RNA (mRNA) molecule during RNA processing to increase the stability of the molecule.
For my project, I would like to estimate the length of poly-A tails from randomly-primed RNAseq data.
Any ideas how this could be done? Are there algorithms/tools which currently do this?
Perhaps I would need to model this somehow? All suggestions appreciated

Why do you have randomly primed RNAseq data? Which organisms are you working with (or it is a mix or a simulated data) ? What have you looked up for your project ? Is there any database of poly-A tail lengths maybe?

@Llopis HeLa cells.  But I would like to do this for a few other organisms as well. "Is there any database of poly-A tail lengths maybe?" It would have to be a data-driven measurement

Just look for polyA tracts at the end of sequences, and count them if they're larger than ~18 bases. I've done this with MinION cDNA reads by mapping the polyA adapter sequence (with an elongated polyA sequence) to the reads using LAST, and working out the length of aligned sequence.
Unfortunately, because the cDNA adapter finishes with a TVN sequence, rather than TTT, binding of the adapter is not necessarily a good representation of the true polyA length. I expect that Illumina adapters will have similar problems, if not worse (due to the read length limits).
It might be possible to use direct-RNA sequencing to better categorise the polyA lengths, but I've yet to look at the public direct-RNA datasets. Just in case you're interested in looking yourself, here's a direct RNA dataset for NA12878 from the wonderful nanopore-WGS consortium:
https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md
Here's the consortium's graph from that summary document showing polyA length:

"> It might be possible to use direct-RNA sequencing to better categorise the polyA lengths, but I'm not familiar with any available public datasets for that." There's no way to infer this from my own dataset?

You've done direct-RNA sequencing on the MinION?

No, I think I'm misunderstanding the reference to "public datasets"

Direct RNA sequencing is a very new technology that can only be done on a MinION. There's not much data available for it yet.

Can’t be done. If you already sequenced then I’m afraid the money is wasted (unless of course the data is good for something else).
The standard Illumina basecaller doesn’t deal well with homopolymers and will consequently call incorrect poly(A) lengths.
Amongst others, Narry Kim’s lab developed a dedicated method — TAIL-seq — to solve this problem. It’s a protocol for library prep using pulldown of biotin tagged tails (rather than random priming), plus dedicated basecalling software that uses spike-ins for calibration.
Be warned though: the protocol is reportedly technically challenging, and in order to use it you need an Illumina machine that gives you access to the raw image files before basecalling. Up-to-date Illumina machines won’t give you these files, full stop. So you’ll need to find a machine that runs an outdated firmware.

(This is a copy of my answer from Reddit; I only now noticed this question here.)

Thanks. I'll probably remove the Reddit question soon, so the link will be broken.


text:
Using Python, how to convert a pandas DataFrame into a VCF?
Let's say I have a pandas dataframe with fields CHROM, POS, ALT, REF. In this special case, I also wouldn't care about IDor FILTER, INFO could be blank (or meaningless) and we'll write QUAL as each 40.
import pandas as pd
example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)
print(df)

   CHROM      POS   REF ALT
0     20    14370     G   A
1     20    17330     T   A
2     20  1110696     A   G
3     20  1230237     T   C
4     20  1234567  GTCT   G

In this case, saving the output with the extension .vcf via df.to_csv("myfile.vcf") with the other required fields would created the tab-delimited format of the VCF, but it wouldn't include the mandatory header. https://samtools.github.io/hts-specs/VCFv4.1.pdf
##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
...
#CHROM POS ID REF ALT QUAL FILTER INFO

The most popular package to work with VCF's appears to be PyVCF: https://pypi.python.org/pypi/PyVCF
However, one needs to begin with a VCF as a template, and add to this (from the documentation):
vcf_reader = vcf.Reader(filename='vcf/test/tb.vcf.gz')
vcf_writer = vcf.Writer(open('/dev/null', 'w'), vcf_reader)
for record in vcf_reader:
     vcf_writer.write_record(record)

It appears the only information necessary from vcf.Reader(filename='vcf/test/tb.vcf.gz') is the metadata?
print(vcf_reader.metadata)
OrderedDict([('fileformat', 'VCFv4.0'), ('fileDate', '20090805'), ('source', ['myImputationProgramV3.1']), ('reference', '1000GenomesPilot-NCBI36'), ('phasing', ['partial'])])

Otherwise, tb.vcf.gz already has variants, which I think means I would need to feed in a VCF with no variants?:
for record in vcf_reader:
     print(record)

Record(CHROM=20, POS=14370, REF=G, ALT=[A])
Record(CHROM=20, POS=17330, REF=T, ALT=[A])
Record(CHROM=20, POS=1110696, REF=A, ALT=[G, T])
Record(CHROM=20, POS=1230237, REF=T, ALT=[None])
Record(CHROM=20, POS=1234567, REF=GTCT, ALT=[G, GTACT])

You will have to define the header either from an existing VCF file or hardcoded into you python script.  Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file.
You will need to order the dataframe columns to match the VCF spec and insert the empty/ default ID, FILTER, INFO and QUAL columns
import pandas as pd

example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)

df['ID'] = ""
df['QUAL'] = 40
df['FILTER'] = ""
df['INFO'] = ""

df = df[['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']]

header = """##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
"""

output_VCF = "myfile.vcf"
with open(output_VCF, 'w') as vcf:
    vcf.write(header)

df.to_csv(output_VCF, sep="\t", mode='a', index=False)

Then myfile.vcf looks like
##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
CHROM   POS ID  REF ALT QUAL    FILTER  INFO
20  14370       G   A   40
20  17330       T   A   40
20  1110696     A   G   40
20  1230237     T   C   40
20  1234567     GTCT    G   40

Thanks for the help. BTW, is fileDate year-month-day or year-day-month? I'm guessing there's no standard.

fileDate is not a required header field you can remove it it you want


text:
What motif finding software is available for multiple sequences ~10Kb?
I have around ~3,000 short sequences of approximately ~10Kb long. What are the best ways to find the motifs among all of these sequences? Is there a certain software/method recommended?
There are several ways to do this. My goal would be to:
(1) Check for motifs repeated within individual sequences
(2) Check for motifs shared among all sequences
(3) Check for the presence of "expected" or known motifs
With respect to #3, I'm also curious if I find e.g. trinucleotide sequences, how does one check the context around these regions?
Thank you for the recommendations/help!

Are you looking for motifs shared by all the sequences? Motifs repeated in each sequence? Are you looking for the presence of specific, known motifs?

"Are you looking for motifs shared by all the sequences? Motifs repeated in each sequence?" I was looking for both motifs shared by all sequences and repeated in sequences, but actually it would be interesting to check each sequence for a repeated sequence. "Are you looking for the presence of specific, known motifs?" No, but it would be interesting to look into this in retrospect

OK, please [edit] your question and ass this information. Comments are easy to miss, hard to read and can be deleted without warning. The three problems (de-novo shared motif identification, de-novo repeated motif identification and detecting known motifs) are different and each requires its own approach.

@terdon Of course

Do you mean 10kb in total, or 3000 sequences, each of length around 10kb?

@CharlesE.Grant 3000 sequences, each of length around 10kb

For (3), this page has a lot of links to pattern/motif finding tools. Following through the YMF link on that page, I came across the University of Washington Motif Discovery section. Of these projection seemed to be the only downloadable tool. I find it interesting how old all these tools are; maybe the introduction of microarrays and NGS has made them all redundant.
Your sub-problem (2) seems similar to the problem I'm having with Nippostrongylus brasiliensis genome sequences, where I'd like to find regions of very high homology (length 500bp to 20kb or more, 95-99% similar) that are repeated throughout the genome. These sequences are killing the assembly.
The main way I can find these regions is by looking at a coverage plot of long nanopore reads mapped to the assembled genome (using GraphMap or BWA). Any regions with substantially higher than median coverage are likely to be shared repeats.
I've played around in the past with chopping up the reads to smaller sizes, which works better for hitting smaller repeated regions that are such a small proportion of most reads that they are never mapped to all the repeated locations. I wrote my own script a while back to chop up reads (for a different purpose), which produces a FASTA/FASTQ file where all reads are exactly the same length. For some unknown reason I took the time to document that script "properly" using POD, so here's a short summary:

Converts all sequences in the input FASTA file to the same length.
     Sequences shorter than the target length are dropped, and sequences longer
     than the target length are split into overlapping subsequences covering
     the entire range. This prepares the sequences for use in an
     overlap-consensus assembler requiring constant-length sequences (such as
     edena).

And here's the syntax:
$ ./normalise_seqlengths.pl -h
Usage:
    ./normalise_seqlengths.pl <reads.fa> [options]

  Options:
    -help
      Only display this help message

    -fraglength
      Target fragment length (in base-pairs, default 2000)

    -overlap
      Minimum overlap length (in base-pairs, default 200)

    -short
      Keep short sequences (shorter than fraglength)

The MEME Suite web site contains a collection of tools for motif analysis (I'm one of the maintainers). It contains two de novo motif discovery tools: MEME and DREME. Public web applications are provided, but you can also download and build command line tools for a local installation.
For your first goal you could use MEME and select the "Any number of repetitions model" (ANR). For your second goal, you'd use MEME with the "Zero or One Occurrences Per Sequence" (ZOOPS) model. For your third goal you could use FIMO (Find Individual Motif Occurrences), and one or more of the motif databases provided on the software and database download page.
It sounds like your sequence data is about 30Mb. The MEME web application is limited to 60kb of sequence data, so you'd have to install a local copy of the MEME Suite. MEME would take a long time to analyze a 30Mb sequence database unless you have MPI configured, and lots of cores available. You might want to consider analyzing a randomly selected subset of your sequences. The running time of MEME grows as the cube of the number of sequences.
For short motifs, you may want to use DREME rather than MEME. DREME is better than MEME at identifying short motifs, but is limited to motifs <= 8 positions wide.

Check out HOMER. "Software for motif discovery and next generation sequencing analysis", it's what my lab uses currently for finding eRNA motifs.
Edit: For @ShanZhengYang "HOMER was designed as a de novo motif discovery algorithm..."
HOMER De Novo Motif

Can one do de novo motif discovery? If not, I'm not entirely sure how I would create a background with this model.

@ShanZhengYang See my edit.

Thanks. I'm still not sure given my set-up how to choose a background for Homer however...

I've never used it for De Novo, so take this with a grain of salt, but try the defaults first. Perhaps shoot the developers an email.

To expand on @EMiller 's comment, for many model organisms, HOMER can figure out a background model for you. You just have to install the supporting files for your model organism using  configureHomer.pl

Most tools I know of looks for enrichment of specific motifs - but that requires that you have a set of sequences which are of special interest and a background set to test against.
Is that your case?
Update after comments 12th June 2017.
You could try the meme suite more specifically the motif finder

Hi Kristoffer Vitting-Seerup, thanks for your comment and welcome to Bioinformatics Stack Exchange. We appreciate any comments like this that help people ask better questions. I have flagged your response as "not an answer" because I think it should be added as a comment to the question (which you can't do yourself until you reach a reputation of 50).

Hi there, yes this would be a better comment to the original question. This is the SO style. Regarding your question, I don't have a background set

it is under development, but maybe BaMMmotif! is something for you? Its main selling point is that it can look for motifs enriched in a set of sequences of equal length de novo. If you can't/don't want to supply a negative set it learns one from the positive sequences. There is a wealth of options to choose from if you have more information about your sequences: there are different models for "zero or one", "one" and "multiple" occurrences of the motif.
You can also use it to look for known motifs, if you encode them as an XXmotif PWM. If you have a file with motifs (like binding sites) you can use that as initialization as well.
While I have not used the software myself, the authors are very responsive on git and the installation instructions seem pretty straightforward.
EDIT: Apparently the software was developed with ChIP experiments as the usual use case scenario, it might misbehave for larger sequences or take super long to run.

Hi galicae, thanks for your answer and welcome to Bioinformatics Stack Exchange. Your answer is a bit short, and could do with a bit more explanation. Could you please add a bit more information into your answer about BaMMmotif? Is it a tool that you have used before? Would this be an appropriate program to find novel shared motifs, repeated motifs within the same sequence, or known motifs?

@gringer fair enough, edited.


text:
Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome?
Let's say I have the genome hg19 loaded into R via BSGenome
library("BSgenome")
hg19genome = getBSgenome('BSgenome.Hsapiens.UCSC.hg19', masked=FALSE)

I then have a list of SNPs loaded as a GRanges object, gr
library(GenomicRanges)

 > gr
 GRanges object with 212 ranges and 3 metadata columns:
        seqnames               ranges strand |     width     REF    ALT
           <Rle>          v <IRanges>  <Rle> | <numeric>
    [1]        1 [86099032, 86099032]      * |         1     C      T
    [2]        1 [86099033, 86099033]      * |         1     C      A
    [3]        1 [86099199, 86099199]      * |         1     G      A
  ....

Is there a straightforward way to inject these SNPs into hg19genome?

By inject you mean that the hg19genome should have the alternative base instead of the reference base for those positions?

@Llopis Yes, that's what I mean

This seems relatively complicated given the structure of a BSGenome object.
The creator of the package answered this question previously on the Bioconductor support forums:
https://support.bioconductor.org/p/86665/#86757

We don't provide an easy way to inject arbitrary SNPs in an arbitrary
  BSgenome at the moment. However, it should not be too hard to forge
  the BSgenome ... package. First you would need to
  compute the sequences of the mutated chromosomes (you can
  use replaceLetterAt for this), then write them to a 2bit file (put
  them in a DNAStringSet object and call rtracklayer::export on it),
  then use that 2bit file to forge the BSgenome ....
  package (see the BSgenomeForge vignette in the BSgenome package for
  how to do this).
  - Herve Pages

Here is another similar question with more info about replaceLetterAt: https://support.bioconductor.org/p/26199/

I noticed those. Question: can you think of another algorithmically efficient way to accomplish this? One way would be to convert the GRanges into a BED, then used this to create a new FASTA with the SNPS, and loaded this into BSGenome and then used getSeq(). That's very complicated though

I should mention that this is a section in the official documentation---the creators appears to have "not gotten around" to implementing this. See " Injecting known SNPs in the chromosome sequences" here: https://bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/GenomeSearching.pdf


text:
How to subset a GRanges via an argument passed into a function?
Let's say I have following example GRanges:
> library(GenomicRanges)
> gr = GRanges(
+     seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
+     ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
+     strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
+     score = 1:10,
+     GC = seq(1, 0, length=10))
> gr
GRanges object with 10 ranges and 2 metadata columns:
    seqnames     ranges strand |     score                GC
       <Rle>  <IRanges>  <Rle> | <integer>         <numeric>
  a     chr1 [101, 111]      - |         1                 1
  b     chr2 [102, 112]      + |         2 0.888888888888889
  c     chr2 [103, 113]      + |         3 0.777777777777778
  d     chr2 [104, 114]      * |         4 0.666666666666667
  e     chr1 [105, 115]      * |         5 0.555555555555556
  f     chr1 [106, 116]      + |         6 0.444444444444444
  g     chr3 [107, 117]      + |         7 0.333333333333333
  h     chr3 [108, 118]      + |         8 0.222222222222222
  i     chr3 [109, 119]      - |         9 0.111111111111111
  j     chr3 [110, 120]      - |        10                 0
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

It is straightforward to subset this GRanges as one would subset an R data.frame, e.g.
> gr$score
 [1]  1  2  3  4  5  6  7  8  9 10

However, I am writing a function where the user would pass through a string for the metadata column to subset e.g.
subsetter = function(granges, column){
    return(granges$column)
}

I try to use the function as follows:
> subsetter(gr, score)

which outputs NULL
This obvious doesn't work, as it will subset the gr by column, not score.
> subsetter(gr, score)
NULL

If this was a data.frame, I could subset as follows:
gr[[(score)]]

but this gives the error
Error in gr[[(score)]] : this S4 class is not subsettable

Are there any other ways to subset a GRanges besides using $?

subsetter = function(gr, cname) {
    return(mcols(gr)[[cname]])
}

You can then use things like subsetter(gr, "GC") and subsetter(gr, "score").

In such cases it is best to check the str of an object gr. Then we could see that meta data is just a dataframe inside S4 object:
gr@elementMetadata
# DataFrame with 10 rows and 2 columns
#        score        GC
#    <integer> <numeric>
# 1          1 1.0000000
# 2          2 0.8888889
# ...

Once we know how to access the dataframe, we can then subset the same way as we would subset a dataframe:
gr@elementMetadata[, "score" ]

I don't see the point of wrapping this into a function, but if it is needed:
subsetter <- function(gr, cname) gr@elementMetadata[, cname ]

# test the function
subsetter(gr, "score")
# [1]  1  2  3  4  5  6  7  8  9 10

Edit
Alternatively, we could use GenomicRanges::elementMetadata() function to extract dataframe then subset:
subsetter <- function(gr, cname) elementMetadata(gr)[, cname ]

While possible, I’d generally avoid accessing S4 object slots: unless explicitly documented, they should be considered an implementation detail that may change without warning. Use the public accessors instead (in this case, mcols).

@KonradRudolph good point, added alternative solution.

Ah … well for that function the documentation says not to use it. ;-)

@KonradRudolph I see, should have read the manuals first. Thank you.

In addition to Devon’s answer, you can use non-standard evaluation to make this work with R symbol names rather than string arguments:
subsetter = function (gr, cname) {
    cname = as.character(substitute(cname))
    mcols(gr)[[cname]]
}

Now using subsetter(gr, score) works analogously to gr$score, which may be desirable. But, likewise, you’ll no longer be able to use a variable:
var = 'score'
subsetter(gr, var) # fails, same as `gr$var`


text:
Nanopore FASTQ header specifications
Is there a specification anywhere for how a Nanopore FASTQ header should be formatted?
I am creating strategies for generating test FASTQ files and would like to ensure I am generating "correct" nanopore sequence identifiers.

Could you clarify? Do you want to simulate ONT reads and make sure the headers are realistic? Or do you have actual reads and want to understand what information the actual headers encode? Or something else?

@DanielStandage see updated question

Thanks. I think that improves the question a lot!

Are you just looking for FASTQ format spec? https://en.wikipedia.org/wiki/FASTQ_format#Format

The headers are defined by the basecaller (guppy, Albacore, bonito, etc) during basecalling. They have changed many times over the years, but generally they start with a readID of @+UUIDv4 followed by multiple entries separated by spaces, of form key=value
For example:
@945da76d-b7ef-4c65-91a1-e69c85832185 runid=2cf9e7c4dd1888d9a1090ebc394dd5b7cb2fe3f0 read=126 ch=181 start_time=2017-10-17T06:43:19Z barcode=barcode01

This is still valid, however, there is some movement to add sam tags to parse data through to sam files when aligned.


text:
Convert FASTA to FASTQ with dummy quality scores
I have a FASTA file which I would like to convert into FASTQ format as the tool I want to use my data in requires it in FASTQ format. So dummy quality scores are fine.
Note: I am not using any functionality of the tool that will require the quality score.

One way of doing this is with two subcommands from the pyfastaq suite.
fastaq to_fake_qual in.fasta - | fastaq fasta_to_fastq in.fasta - out.fastq

The first tool, to_fake_qual, creates fake quality scores (default 40) for each base and the - sends that file (.qual) to stdout. The second tool, fasta_to_fastq, consumes both the original fasta and the quality scores coming from stdin and turns these into a fastq file.

A simple Biopython solution:
from Bio import SeqIO

for r in SeqIO.parse("myfile.fa", "fasta"):
    r.letter_annotations["solexa_quality"] = [40] * len(r)
    print(r.format("fastq"), end='')

Example:
$ cat << EOF > myfile.fa
>1
ATG
>2
ATGTAGA
EOF
$ python3 myscript.py
@1
ATG
+
III
@2
ATGTAGA
+
IIIIIII

Here's a quick one-liner with no dependencies other than the UNIX core utils and Perl.
cat seqs.fasta | paste - - | perl -ne 'chomp; s/^>/@/; @v = split /\t/; printf("%s\n%s\n+\n%s\n", $v[0], $v[1], "B"x length($v[1]))' > seqs.fastq

Caveat: assumes your Fasta sequences are not wrapped, i.e. two lines per record.

I think the python package bioconvert does what you want if you don't provide qualities: https://github.com/bioconvert/bioconvert/blob/master/bioconvert/fasta2fastq.py
The following outputs lots of warnings for me, but does generate a fastq file test.fq from an existing fasta file test.fa (using dummy quality "I", and adding "None" to sequence headers that do not have a comment part in the original fasta):
bioconvert fasta2fastq test.fa test.fq

(This requires installing bioconvert using pip beforehand: pip install bioconvert.)

I have just uploaded a Python tool to convert fasta to fastq: https://github.com/Mostafa-MR/A_tool_to_convert_fasta_to_fastq
Please note this is presented directly from a Jupyter notebook, i.e. it is a ipynb.


text:
How to convert fastq to fast5
fast5 is a variant of HDF5 the native format in which raw data from Oxford Nanopore MinION are provided. You can easily extract the reads in fast5 format into a standard fastq format, using for example poretools.
Say I have aligned these reads in fastq format to an external reference genome, resulting in a SAM file. Say I have then taken a subset of the  SAM file, according to the bitwise flag, to include only the reads that map to the reference. With the read ID, I can then grep them out from the file containing the reads in fastq format, generating a subset file in fastq format containing only the IDs that have mapped to the reference.
Now my question is, can we subset reads from the fast5 archive according to the list of mapping reads as taken from the file with reads in fastq format? This is for educational purposes, so that we have a smaller starting archive, and the fast5 -> fastq extraction takes less cpu time.

This seems a nice question to play around if we could have sample data. Could you provide (or a link to) a fast5 and the list of mappings ?

There are some sample reads from my TEDx talk last year here, reads from runs I did during the AGTA conference last year here, and reads from a chimeric read investigation here. The ONT FAST5 "standard" change so frequently that it's difficult to find something that represents what's currently available.

NOTICE:
I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same.

This is something I have been meaning to get around to for a while, so thanks for the prompt.
I have created a python program called fast5seek to do what (I think) you're after.
As you have mentioned this is for educational purposes I have added a tonne of comments to the code too so I think you shouldn't have any issues following it.
The docs on the GitHub repo have all the info, but for those reading along at home
pip3 install fast5seek
fast5seek -i /path/to/fast5s -r in.fastq in.bam in.sam -o out.txt

What it does is read in <in.fastq|in.bam|in.sam> and extract the read id from each header. It then goes through all the fast5 files under /path/to/fast5s and checks whether their read id is in the set of read ids from <in.fastq|in.bam|in.sam>. If it is, the path to the file is written to it's own line in out.txt.
If no output (-o) is given, it will write the output to stdout.
So if you wanted to pipe these paths into another program, you could do something like
mkdir subset_dir/
fast5seek -i /path/to/fast5s/ -r in.fastq | xargs cp -t subset_dir/

The above example would copy the fast5 files that are found in your fastq/BAM/SAM to subset_dir/.

As a caution, the run id should be checked (if that information is available). There's a small possibility that two FASTQ files from different runs will have the same read id.

I have now added this check for fastq files. Although, run_id is not contained within BAM/SAM files (at least not for the ones I have) so this check will only occur on fastq.

I do this very frequently, using the read file name to identify FAST5 files associated with particular reads. If a FASTQ record includes the channel number and read number (and preferably the runID as well), then I use that information to find the associated FAST5 read.
If reads are called with Albacore, then the sequencing_summary.txt file has additional information. The first column is the exact name of the file that was called, the second column is the read ID of the FASTQ sequence, and the third column is the run ID associated with that read ID. This is much easier to work with, but does require calling reads with Albacore (which, admittedly, is what seems to be producing the best results at the moment).

yes this is what I did, basically: 1. extract fast5 to fastq; 2. align the fastq to the reference; 3. get the reads ID that map from the SAM of alignment; 4. grep the mapping reads ID from the fastq; 5. extract the name of the fast5 from the fasta header of each read (last field); 6. find the corresponding fast5 in the original fast5 directory. But I was wondering if I did not miss any software that could do this on one run.

Added a note about albacore, which has a direct link between FASTQ read IDs and FAST5 file names

In light of your above comment @AminaEchchiki I have added in compatibility for BAM and SAM files too.


text:
simple data subsetting on r
I have a big df like following (pseudo) df and I want to retrieve values in different col. according to the rows in the first col.
A B C D
a 1 1 1
a 2 3 1
a 3 3 3
a 4 5 6
b 5 6 5
b 7 6 5
c 7 6 5
c 5 4 3
c 5 5 4

I want to retrieve the values in col. B, C, and D according to the rows in col A and perform further analysis of them. For example, multiply the values in col B and C for all the rows of "a" in col A. Can someone suggest a commend for such a naive question. Thanks.

I'm downvoting this because it doesn't have anything to do with bioinformatics, but you can use df[df$A=='a',], or dplyr::group_by(A) if you want to repeat the same analysis for each value of A

That is not true, the original df contains genetic info. Thanks

I'm voting to close this question as off-topic because without more context, this looks like a programming question better suited for Stack Overflow, and the "this question belongs on another site in the SE network" radio button does not offer Stack Overflow as a selection.

Here's an approach that uses native R, using the vectorized (fast) family of apply() functions.
> df2 <- df[df$A == 'a', which(names(df) %in% c("B","C"))]
> apply(df2, 1, prod)

Use apply/lapply/mapply/etc. where you can in place of alternatives, where you don't know their run characteristics.
Here's a more concrete example:
> df
  A B C D
1 a 1 1 1
2 a 2 3 1
3 a 3 3 3
4 a 4 5 6
5 b 5 6 5
...
> df2 <- df[df$A == 'a', which(names(df) %in% c("B","C"))]
> df2
  B C
1 1 1
2 2 3
3 3 3
4 4 5
> apply(df2, 1, prod)
 1  2  3  4
 1  6  9 20

Edit: If you need to do this on each unique value of df$A, you can nest apply calls to make this scale:
> perA.list <- apply(data.frame(unique(df$A)), 1, function(x) apply(df[df$A == x, which(names(df) %in% c("B","C"))], 1, prod))
> perA.list
[[1]]
 1  2  3  4
 1  6  9 20

[[2]]
 5  6
30 35

[[3]]
 7  8  9
35 20 25

Each index in this list correspondents to the index of the value from unique(df$A), i.e. a, b, and c.
If you don't want to do this for all unique values in df$A, but for some specific element in df$A, then you would just need to write a function. 95% of the code is in this answer to write such a function, so I'll leave that as a homework exercise for you to implement.

Thanks, I can't able to mention "a"every time. First I also want to have same results for "b" and "c" as for "a". Second, as I mentioned the data is big with 70,000 entries. In this case please suggest something.

See edited answer.

Hi, ddply (https://www.rdocumentation.org/packages/plyr/versions/1.8.4/topics/ddply) was used to obtain results, exactly what I was looking. Thanks


text:
sum of products of mapq and mapped bases for each read in a from a BAM file
Given a BAM file, I'd like to calculate the sum, over all reads, of the mapping quality and the number of mapped bases (i.e. number of M's in the CIGAR string).
For example, given two reads like this:
read1: mapq=40, mapped_bases=10
read2: mapq=20, mapped_bases=50

The number I'm looking for is:
(40*10) + (20*50) = 1400

One option here is to stream the file and get the two quantities line-by-line. But this seems likely to be slow. Thoughts welcome.

BAM files are stored line by line, so your only option is to process it that way. In python with pysam:
import pysam

bam = pysam.AlignmentFile("something.bam")
tot = 0
for read in bam.fetch():
    if read.is_unmapped:
        continue
    tl = 0
    for op, l in read.cigartuples:
        if op in [0, 7, 8]:
            tl += l
    tot += tl * read.mapping_quality
bam.close()
print("the total is {}".format(tot))

One possibility uses bamutils (http://ngsutils.org/modules/bamutils/removeclipping/) to remove the clipped bases, then awk to get the sum of the products, something like:
bamutils removeclipping val.bam val_rc.bam
samtools view val_rc.bam | awk '{print length($10)*$5}' | awk '{sum+=$0} END { print "mapq*len sum =",sum}'

the bamutils step is rather slow, but it will do the job.

Why use three separate awk calls instead of samtools view val_rc.bam | awk '{sum+=length($10)*$5} END { print "mapq*len sum =",sum}' ?


text:
Minion channel ID's from Albacore
The latest version of Albacore from Oxford Nanopore Technologies calls bases from raw fast5 files. A useful piece of output is the sequence_summary.txt, which is a big tab-delimited file with information on each read.
One of the columns in there is the channel, which references the channel number from 1 to 512 on the flowcell.
These channels obviously have a physical layout on the flowcell. What is the layout? I.e. how many channels are there per row and per column, and how do the channel numbers map to the rows and columns?
My best guess is that it's a 16*32 array, but I don't know for sure, and even if I did I wouldn't know whether the channels map across the shorter or the longer axis in order (or maybe they are not in order at all, but that would be silly).

thanks to the selected answer, we have it figured out. In case anyone needs an implementation in R, here's one: https://gist.github.com/roblanf/df47b9748c3aae00809cc675aca79989

There are actually 2048 usable sequencing wells, hexagonally packed with four wells connected to the same sequencing sensor/channel via a multiplex (mux) selector. The combination of the mux and the channel number determines the physical location of the well. Unfortunately, the association between channel number and physical location is not obvious, and within each channel the muxes are also not in an obvious order, i.e. [3,4,1,2,2,1,4,3] for two adjacent channels.
Here's an "official" text description of the layout as it was for R7; it hasn't changed for the most recent flow cells (R9.5):

Channels 1-64 occur at the top of the chip (the other low channel numbers are at the bottom)
Channels order down the chip: 1-64, 449-512, 385-448, 321-384, 257-320, 193-256, 129-192, 65-128
Muxes run from left to right in the order: 3, 4, 1, 2, 2, 1, 4, 3

To confirm this layout, I pulled the channel/position lookup table definitions out of MinKNOW, and it seems to be the same as this:

Here's the code to produce this image (I haven't adjusted it to have a hexagonal structure):
#!/usr/bin/Rscript
library(rjson);

layout.fname <- "/opt/ONT/MinKNOW/Client/resources/app/resources/channels.json";

data.list <- fromJSON(file=layout.fname);

png("MinION_FC_Layout.png", width=3072, height=1536, pointsize=24);
par(mar=c(4,4,0.5,0.5), lwd=3);
data.mat <- matrix(unlist(data.list), ncol=2, byrow=TRUE);
plot(data.mat, pch=21,
     bg=(colorRampPalette(c("white","grey"))(2048))[1:2048],
     col=c("red","yellow","green","blue"),
     xlab="locX", ylab="locY", cex=2,
     xlim=c(7,76), ylim=c(0,33));
text(x=data.mat, labels=rep(1:512, each=4), cex=0.4);
legend("top",legend=c(1,2,3,4), fill=c("red","yellow","green","blue"),
       inset=0.05, title="MUX");
text(x=6, y=data.mat[data.mat[,1] == 8,2],
     labels = ceiling(which(data.mat[,1] == 8) / 4));
text(x=77, y=data.mat[data.mat[,1] == 75,2],
     labels = ceiling(which(data.mat[,1] == 75) / 4));
arrows(x0=c(8,75), x1=c(39,44), y0=33, length=0.2);
invisible(dev.off());

an interesting thing to note is that this answer is different from @wouter-de-coster answer that comes from poRe. As far as I can tell, this is the correct answer.

I think the layout in poRe was based on the numbering in the original MinKNOW GUI (which was also incorrect).

for those interested I built an annoying by-hand implementation of the mapping to use for plotting in R here: https://gist.github.com/roblanf/df47b9748c3aae00809cc675aca79989

poRe has a show.layout() function which shows you the 32*16 grid on which the channels are arranged. EDIT: THIS APPEARS NOT TO BE THE CORRECT PHYSICAL LAYOUT.
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
 [1,]  125  121  117  113  109  105  101   97   93    89    85    81    77
 [2,]  126  122  118  114  110  106  102   98   94    90    86    82    78
 [3,]  127  123  119  115  111  107  103   99   95    91    87    83    79
 [4,]  128  124  120  116  112  108  104  100   96    92    88    84    80
 [5,]  253  249  245  241  237  233  229  225  221   217   213   209   205
 [6,]  254  250  246  242  238  234  230  226  222   218   214   210   206
 [7,]  255  251  247  243  239  235  231  227  223   219   215   211   207
 [8,]  256  252  248  244  240  236  232  228  224   220   216   212   208
 [9,]  381  377  373  369  365  361  357  353  349   345   341   337   333
[10,]  382  378  374  370  366  362  358  354  350   346   342   338   334
[11,]  383  379  375  371  367  363  359  355  351   347   343   339   335
[12,]  384  380  376  372  368  364  360  356  352   348   344   340   336
[13,]  509  505  501  497  493  489  485  481  477   473   469   465   461
[14,]  510  506  502  498  494  490  486  482  478   474   470   466   462
[15,]  511  507  503  499  495  491  487  483  479   475   471   467   463
[16,]  512  508  504  500  496  492  488  484  480   476   472   468   464
      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
 [1,]    73    69    65    61    57    53    49    45    41    37    33    29
 [2,]    74    70    66    62    58    54    50    46    42    38    34    30
 [3,]    75    71    67    63    59    55    51    47    43    39    35    31
 [4,]    76    72    68    64    60    56    52    48    44    40    36    32
 [5,]   201   197   193   189   185   181   177   173   169   165   161   157
 [6,]   202   198   194   190   186   182   178   174   170   166   162   158
 [7,]   203   199   195   191   187   183   179   175   171   167   163   159
 [8,]   204   200   196   192   188   184   180   176   172   168   164   160
 [9,]   329   325   321   317   313   309   305   301   297   293   289   285
[10,]   330   326   322   318   314   310   306   302   298   294   290   286
[11,]   331   327   323   319   315   311   307   303   299   295   291   287
[12,]   332   328   324   320   316   312   308   304   300   296   292   288
[13,]   457   453   449   445   441   437   433   429   425   421   417   413
[14,]   458   454   450   446   442   438   434   430   426   422   418   414
[15,]   459   455   451   447   443   439   435   431   427   423   419   415
[16,]   460   456   452   448   444   440   436   432   428   424   420   416
      [,26] [,27] [,28] [,29] [,30] [,31] [,32]
 [1,]    25    21    17    13     9     5     1
 [2,]    26    22    18    14    10     6     2
 [3,]    27    23    19    15    11     7     3
 [4,]    28    24    20    16    12     8     4
 [5,]   153   149   145   141   137   133   129
 [6,]   154   150   146   142   138   134   130
 [7,]   155   151   147   143   139   135   131
 [8,]   156   152   148   144   140   136   132
 [9,]   281   277   273   269   265   261   257
[10,]   282   278   274   270   266   262   258
[11,]   283   279   275   271   267   263   259
[12,]   284   280   276   272   268   264   260
[13,]   409   405   401   397   393   389   385
[14,]   410   406   402   398   394   390   386
[15,]   411   407   403   399   395   391   387
[16,]   412   408   404   400   396   392   388

this solution (while I can see that it should be thought to be correct) appears to be at odds with the current information (in the selected answer). I wonder if poRe uses a layout from an older flowcell?

Hmm good catch.

for those interested, I wrote a quite compact Python function to generate the (correct!) layout in a numpy array: https://gist.github.com/wdecoster/1a86d4ba1039afebb250352fc4e502d1


text:
Confirm success or failure of RNA-Seq normalization
I am working with a set of (bulk) RNA-Seq data collected across multiple runs, run at different times of the year. I have normalized my data using library size / quantile / RUV normalization, and would like to check (quantitatively and/or qualitatively) whether or not normalization has succeeded in removing the batch effects.
It is important to note that by "normalization has succeeded", I simply mean that unwanted variation has been removed - further analysis is required to ensure that biological variation has not been removed. What are some plots / statistical tests / software packages which provide a first-step QC for normalization?

You should use box plots and PCA plot. Let's take a look at the RUV paper:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4404308/

Before normalization and after UQ normalization:

Libraries do not cluster as expected according to treatment. ... for UQ-normalized counts. UQ normalization does not lead to better clustering of the samples...

Before normalization, the medians in the box-plot obviously look very different among replicates.
After UQ normalization, the medians look closer but Trt.11 look like an outlier. Furthermore, the treatments aren't clustered on the PCA plot. Since they are replicates, you'd like them be close on the plot.
After RUV normalization

... RUVg shrinks the expression measures for Library 11 towards the median across libraries, suggesting robustness against outliers. ... Libraries cluster as expected by treatment. ...

The RUV has made the distribution more robust and the samples closer on the PCA plot. However, it's still not perfect as one of the treatments is not close to the other two on the first PC.
The vignettes for Bioconductor RUVSeq describes the two functions: plotRLE and plotPCA.

I might point out the distinction between an RLE plot (shown here) and an ordinary boxplot (a distinction made on page 3 of the April 2017 preprint on RLE Plots: Visualizing Unwanted Variation in High Dimensional Data at https://arxiv.org/pdf/1704.03590.pdf). Otherwise an excellent answer but important to explain what RLE means.

Visual inspection with histograms, boxplots, or some other distribution visualization is the way to go. Prior to normalization, your abundances may look something like this.

Post-normalization, they should look something like this.

See this blog post for example code.


text:
Uppercase vs lowercase letters in reference genome
I am using a reference genome for mm10 mouse downloaded from NCBI, and would like to understand in greater detail the difference between lowercase and uppercase letters, which make up roughly equal parts of the genome. I understand that N is used for 'hard masking' (areas in the genome that could not be assembled) and lowercase letters for 'soft masking' in repeat regions.

What does this soft masking actually mean?
How confident can I be about the sequence in these regions?
What does a lowercase n represent?

What does this soft masking actually mean?

A lot of the sequence in genomes are repetitive. Human genome, for example, has (at least) two-third repetitive elements.[1].
These repetitive elements are soft-masked by converting the upper case letters to lower case. An important use-case of these soft-masked bases will be in homology searches: An atatatatatat will tend to appear both in human and mouse genomes but is likely non-homologous.

How confident can I be about the sequence in these regions?

As you can be about in non soft-masked based positions. Soft-masking is done after determining portions in the genome that are likely repetitive. There is no uncertainty whether a particular base is 'A' or 'G', just that it is part of a repeat and hence should be represented as an 'a'.

What does a lowercase n represent?

UCSC uses Tandom Repeat Finder and RepeatMasker for soft-masking potential repeats. NCBI most likely uses TANTAN. 'N's represents no sequence information is available for that base. It being replaced by 'n' is likely an artifact of the repeat-masking software where it soft-masks an 'N' by an 'n' to indicate that portion of the genome is likely a repeat too.
[1] http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002384

Informative answer, but I think it's controversial to say the human genome is "(at least) two-third repetitive elements"; the P-clouds method you cite is quite permissive and half is a more commonly accepted figure. And soft-masking doesn't involve masking all repeats generally, just interspersed repeats and low complexity sequences. Also there is always uncertainty around base calling and assembly building, and more so for repetitive sequences, although mm10 is one of the best assemblies of course...

The use of lower/upper case letters and N/n letters in genomes sequences is not completely standardised and you should always check the specification of the resource you are using.
Lower case letters are most commonly used to represent “soft-masked sequences”, a convention popularised by RepeatMasker, where interspersed repeats (which covers transposons, retrotransposons and processed pseudogenes) and low complexity sequences are marked with lower case letters. Note that larger repeats, such as sizable tandem repeats, segmental duplications, and whole gene duplications are not generally masked.
However, there are other uses for lower/upper case letters, for example, Ensembl have used upper/lower case letters to represent exonic and intronic sequences respectively.
N and n nucleotides may represent “hard masked sequences”, where interspersed repeats and low complexity sequences are replaced by Ns. But N/ns may alternatively represent ambiguous nucleotides, indeed this is the IUPAC specification.
Also note occasionally (although fortunately rarely) X/x is used to represent ambiguous nucleotides or “hard-masked sequences” too.

Lower case nucleotides commonly denotes a soft masked sequences. How exactly the genome was masked you can find in FAQ of NCBI:

Are repetitive sequences in eukaryotic genomes masked?
Repetitive sequences in eukaryotic genome assembly sequence files, as
  identified by WindowMasker, have been masked to lower-case.
The location and identity of repeats found by RepeatMasker are also
  provided in a separate file. These spans could be used to mask the
  genomic sequences if desired. Be aware, however, that many less
  studied organisms do not have good repeat libraries available for
  RepeatMasker to use.

IMHO, low complexity regions are always more likely to be missassembled than high complexity sequences. However, this will be problem for non-model organisms. I would guess that the reliability of the softmasked regions of mouse genome will very high.
No idea, looks like an artefact.

An example of usage of the soft mask
Mapping of sequence to reference usually starts with perfect matches of seeds (substrings) of the mapped reads and the reference sequence. Soft masked (low complexity) regions are not used for matches of seeds, but they are used only for the extension of the alignment if there was a seed in a neighbouring region. This application of softmasking applied to problem of long read assembly is described on this blog.


text:
How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq?
In single-cell RNA-seq data we have an inflated number of 0 (or near-zero) counts due to low mRNA capture rate and other inefficiencies.
How can we decide which genes are 0 due to gene dropout (lack of measurement sensitivity), and which are genuinely not expressed in the cell?
Deeper sequencing does not solve this problem as shown on the below saturation curve of 10x Chromium data:

Also see Hicks et al. (2017) for a discussion of the problem:

Zero can arise in two ways:

the gene was not expressing any RNA (referred to as structural zeros) or
the RNA in the cell was not detected due to limitations of current experimental protocols (referred to as dropouts)

Possible duplicate of What methods are available to find a cutoff value for non-expressed genes in RNA-seq?

I don't think this is a duplicate: this question is about distinguishing between zero counts and the linked question is about distinguishing between zero and non-zero counts.

Definitely not a duplicate, this is a very different question.

Sorry, I missed the "single cell" mention. Despite that, I expect that some answers are likely to overlap, even if the problem is different. See, for example, Kristoffer's answer.... er, which may have been posted in the wrong question.

hey @Peter if the answers here are not satisfactory you could maybe add more info to the question or comment on the answers

Actually this is one of the main problems you have when analyzing scRNA-seq data, and there is no established method for dealing with this. Different (dedicated) algorithms deal with it in different ways, but mostly you rely on how good the error modelling of your software is (a great read is the review by Wagner, Regev & Yosef, esp. the section on "False negatives and overamplification"). There are a couple of options:

You can impute values, i.e. fill in the gaps on technical zeros. CIDR and scImpute do it directly. MAGIC and ZIFA project cells into a lower-dimensional space and use their similarity there to decide how to fill in the blanks.
Some people straight up exclude genes that are expressed in very low numbers. I can't give you citations off the top of my head, but many trajectory inference algorithms like monocle2 and SLICER have heuristics to choose informative genes for their analysis.
If the method you use for analysis doesn't model gene expression explicitly but uses some other distance method to quantify similarity between cells (like cosine distance, euclidean distance, correlation), then the noise introduced by dropout can be covered by the signal of genes that are highly expressed. Note that this is dangerous, as genes that are highly expressed are not necessarily informative.
ERCC spike ins can help you reduce technical noise, but I am not familiar with the Chromium protocol so maybe it doesn't apply there (?)

since we are speaking about noise, you might consider using a protocol with unique molecular identifiers. They remove the amplification errors almost completely, at least for the transcripts that you capture...
EDIT: Also, I would highly recommend using something more advanced than PCA to do the analysis. Software like the above-mentioned Monocle or destiny is easy to operate and increases the power of your analysis considerably.

Some people use imputation to differentiate between true zeros and dropout in single-cell data. Some approaches you can look into:

CIDR which was published in Genome Biology, here
scImpute , published in Nature, here
Seurat (addImputedScore) , i.e. the famous R package, also here
SAVER which was reviewed in Frontiers of Genetics, here