# Question
Title: nextflow: tumor normal sample - how to make code organic
I'm developing pipeline in nextflow it uses tools: fastp, bwa mem index, bwamem, gatk mark duplicates, gatk setupnmd, gatk applybqsr, gatk recalibrate. I do this using FASTQ.gz of normal and tumor DNA samples.
I'm interested to use mutect2, lancet and strelka in nextsteps. However, I am unable to write code in the pipeline that makes sense. I use loads of shell scripting in mutect2 to use tumor-normal.
To expand, I've input in FASTQ as tiny\_n, tiny\_t, where \_n is normal and \_t is tumor. I run pipeline until applybqsr for each individually. However, as I reach mutect2, I check if sample contains **\_n** or otherwise, I then create an empty file for the output as required in the mutect2 process.

How do I overcome creating dummy/empty file?

Please see below code of **main.nf**:

```
params.raw_files="/sc/arion/projects/tiny/all_tumors_normals/tiny_*{R1,R2}.fastq"

include { FASTP} from './fastp_process.nf'
//include {bwa_index} from './index_process.nf'
include { align_bwa_mem} from './bwamem_process_already_index.nf'
include { gatk_markduplicates} from './gatk_markduplicates_process.nf'
include { recalibrate_bam } from './recalibratebam_process.nf'
include {setupnmdtags} from './setupnmdtags_process.nf'
include { applybqsr } from './applybqsr_process.nf'
include { mutect2 } from './mutect2_process.nf'
workflow {

        read_pairs_ch = Channel.fromFilePairs(params.raw_files)

        FASTP(read_pairs_ch)

        align_bwa_mem(FASTP.out.reads)
        gatk_markduplicates(align_bwa_mem.out.sorted_bams)
        setupnmdtags(gatk_markduplicates.out.sorted_bam_mark_duplicates)
        recalibrate_bam(setupnmdtags.out.setupnmdtags_bam)

applybqsr(gatk_markduplicates.out.sorted_bam_mark_duplicates,recalibrate_bam.out.recalib_table)

mutect2(applybqsr.out.recal_bam)
}

```

below is code of mutect2 to handle only one sample, run for normal, create dummy/empty file for tumor sample, **mutect2\_process.nf**:

```
process mutect2 {
        maxForks 3
        debug true
        tag "${sample_id}"
 memory '15 GB'

        input:
        tuple val(sample_id),path(bqsrbam)

        output:
        tuple val(sample_id),path("${sample_id}.mutect2out.vcf")

        script:
                         def avail_mem = task.memory ? task.memory.toGiga() : 0
                def java_options = [
        avail_mem ? "-Xmx${avail_mem}G" : "",
        "-XX:+UseSerialGC",
    ]

        """
        echo "${sample_id}\n"

        grep_value=\\\$(echo "${sample_id}|  | grep "_n" ")

        if [[ ! -z "\\\${grep_value}" ]]
        then

                identifier=\\\$(echo ${sample_id} | awk -F"_" '{print \\\$1}')
                echo "\\\${identifier}"_t.sorted.markdup.recal.bam""

                ml gatk/4.1.3.0
gatk --java-options "${java_options.join(' ')}" Mutect2 -L $params.target -I "${params.outdir}/applybqsr/""\\\${identifier}"_t.sorted.markdup.recal.bam"" -I  "${params.outdir}/applybqsr/""\\\${identifier}"_n.sorted.markdup.recal.bam"" -normal \\\${identifier}"_n" -tumor \\\${identifier}"_t" --output ${sample_id}".mutect2out.vcf" -R $params.hg38genome

        else
                echo "Nothing to do\n"
                touch ${sample_id}".mutect2out.vcf"

        fi

        """
        }

workflow.onComplete {
        log.info ( workflow.success ? "Done mutect2!" : "Oops .. something went wrong in mutect2" )
}

```

There is a load of unix/shell code that does trivial checks. It runs gatk mutect2 oly when \_n is found. How do I overcome: output defined `tuple val(sample_id),path("${sample_id}.mutect2out.vcf")` ?
Is there a better way to do this?

# Answer
I think you have a design problem here. A better way would be perform this operation in your *workflow* block. This would allow you to feed in the tumor and normal BAMs (and avoid accessing files outside of the working directory as previously discussed). You can use the `map` operator to extract the patient or case id from the sample names and the `filter` operator to separate the samples into two channels. You could also use the `branch` operator for this like in the example below. The `join` operator can then be used to join together the items emitted in a consistent way. For example:

```
workflow {

    ref_fasta = file( params.ref_fasta )
    targets = file( params.targets )

    ...

    applybqsr.out.recal_bam
        .branch { sample_id, bam ->

            def case_id = sample_id.substring(0, sample_id.lastIndexOf('_'))

            tumor: sample_id.toUpperCase().endsWith('_T')
                return tuple( case_id, sample_id, bam )

            normal: sample_id.toUpperCase().endsWith('_N')
                return tuple( case_id, sample_id, bam )

        }
        .set { cases_ch }

    joined_ch = cases_ch.tumor.join( cases_ch.normal )

    mutect2( joined_ch, ref_fasta, targets )
}

```

```
process mutect2 {

    tag { sample_id }

    module 'gatk/4.1.3.0'

    input:
    tuple val(case_id), val(tumor), path(tumor_bam), val(normal), path(normal_bam)
    path reference_fasta
    path intervals_bed

    output:
    tuple val(case_id), path("${case_id}.mutect2.vcf")

    script:
    def avail_mem = task.memory ? task.memory.toGiga() : 0
    def java_options = [
        avail_mem ? "-Xmx${avail_mem}G" : "",
        "-XX:+UseSerialGC",
    ]

    """
    gatk \\
        --java-options "${java_options.join(' ')}" \\
    Mutect2 \\
        -R "${reference_fasta}" \\
        -L "${intervals_bed}" \\
        -I "${tumor_bam}" \\
        -I "${normal_bam}" \\
        -normal "${normal}" \\
        -tumor "${tumor}" \\
        -O "${case_id}".mutect2.vcf"
    """
}

```

There's also likely an issue with how *applybqsr* is invoked. If you require multiple input channels, most of the time what you'll want is one queue channel and one or more value channels. I suspect the inputs for *applybqsr* also need to be joined into a single queue channel.
> 2 votes

---
Tags: nextflow
---


# Question
Title: Finding linkage disequilibrium when I have a kinship matrix but not pedigree
I have a bunch of samples of a non-model organism divided in two groups: parental generation and offspring. I do not have a pedigree so I do not actually know who is a parent of who, and which are the sibships. I have a bunch of SNPs that I can use to create a kinship matrix, but I cannot resolve families just using it.

I want to perform linkage disequilibrium (LD) between the SNPs. I can calculate the correlation between them, but I feel I am losing a lot of power because there is kinship between samples. Is there a program that could take the kinship matrix in account when calculating the LD?

I could, in theory, reconstruct the pedigree from the markers. I run multiple instance of Colony but the results are underwhelming -- most samples are not allocated (admittedly this raises questions on the SNP situation, but this is a whole different issue).

# Answer
> 2 votes
I originally thought you were working on humans. I spotted "non model organism", that changes things alot, albeit I guess humans are a non-model organism in a way. The question is what the organism. LD has been applied to every kingdom of life, except viruses. Arlequin isn't bad for relatively small to medium data sets (but tricky to get the format right - thats a separate question). If it's a genome sequence then I still suspect its plink. There are a number of recent LD programs BTW.

I did notice Colony was written in Fortran! Wow, the author must be a physicist at ZSL anyway...

If you have a defined number of SNPs - sounds small to medium data set - then I'd use Arlquein http://cmpg.unibe.ch/software/arlequin35/Arl35Methods.html, Arlquin does loads of pop gen BTW its very famous for AMOVA BTW (very cool calculation). HAPLOVIEW has been used for LD, these authors used it here, what surprised me was they avoided using Arlequin for the LD but used HAPLOVIEW. There's a plethora of recent LD algorithms, never used them so can' say, e.g. quickLD, stuff in R.

*What I am definitely saying is the LD works on from raw alleles* Thats the most important part of this response. I've attempted to describe the rationale below.

---

If this is human genomes, or large genome data I'd would recommend plink analysis performed directly on your data, as described here https://www.cog-genomics.org/plink/1.9/ld

Using any sort of matrix will remove the linkage disequilibrium signal, **IF** its a single matrix for all your data. If its multiple matrices you can perform a correlation statistic between the matrices - but its better just to do a straight linkage disequilibrium calculation, e.g. plink.

The reason why correlation matrices will fail is because LD needs to identify incongruence/disparity and whether that incongruence occurs within random expectations.

Basically you are making a phylogenetic tree - thats what the correlation matrix is about, but the tree assumes clonality. LD is about deviation from random expectations (equilibrium), so either a random distribution, e.g. Poison or a randomisation process is needed. What I'm attempting to say is "a tree" (single matrix) and LD are very different calculations. Whilst, technically you can do LD via matrices (plural), its better to do a LD calculation from the raw alleles.

---

From the comments, it appears a "tree like structure" (ancestors) is wanted to direct an LD calculation.

What you can do is subdivide the populations and perform LD on each subdivision. A tree (here using a correlation or kinship matrix) is a continuum so it's at what point do you call the the subdivision? Its doable but the information is not presented here.

What you can do is use "Structure" - this has been upgraded alot. This is a HWE, I think using Bayes and it designates how many populations (at LD) are present in your data. You can then perform a specific LD directly on each subdivision.

If you compare the output of "Structure",basically Bayesian HWE you should see a correlation with your kinship matrix. What you are inferring is a partial LD parameter that will define ancestry and whilst its doable it has never been calculated in this manner that I'm aware of: its random or it ain't and thats the criteria of significance. This is definitely how LD works in small genomes (everything below higher eukaryotes).

---
Tags: phylogenetics, phylogeny, plink, pedigree
---


title: Help with error DiffBind

text:
# Question
Title: Help with error DiffBind
*This question was also asked on Biostars and the Bioconductor Forum*

I want to find differential accessibility between conditions (WS, WF, MS, MF). I can compare any pair conditions without error except this pair MS vs MF. Data is the bam and broadpeak files from ATAC seq pipeline and the code is like a pipeline of the DiffBind R package. The sample is human and uses Illumina sequencing technology.

Would you please suggest what wrong in this case? I pretty much run the same code as I used for another data set and didn't have this error. Thank you so much!

```
result <- dba(sampleSheet=samples)
result <- dba.blacklist(result)
result <- dba.count(result, bParallel=FALSE)
result <- dba.count(result)
result <- dba.normalize(result)
result <- dba.contrast(result, minMembers=2)
result <- dba.analyze(result)
result <- dba.plotProfile(result)

Generating report-based DBA object... Error: No valid contrasts/methods specified.

```

There is post about this error but I still got error when trying those solutions. Such as the order of the function, the path to the bam and broad peak files in sample sheet.

```
dba.show(result, bContrasts=TRUE)
     Factor    Group Samples  Group2 Samples2 DB.DESeq2
1 Condition diseased       2 control        2         0

```

I don't know why the value of DB.DESeq2 in this case was 0 which maybe the reason of the error.

This is the result after running dba.contrast:

```
4 Samples, 89585 sites in matrix:
   ID Tissue   Factor Condition Replicate    Reads FRiP
1  MS    non ATAC-seq   control         1 23731398 0.17
2 MSR    non ATAC-seq   control         2 10499000 0.16
3  MF    non ATAC-seq  diseased         1 31535510 0.18
4 MFR    non ATAC-seq  diseased         2 18488559 0.18

Design: [~Condition] | 1 Contrast:
     Factor    Group Samples  Group2 Samples2
1 Condition diseased       2 control        2

```

If I try to compare MS vs MF, I will get that error but if I remove MS and replace by other sample, I don't have that error anymore.

I don't have that error if I add more samples and the value of DB.SESeq2 at the first row from dba.show() is not 0.

Full samplesheet:

```
SampleID,Tissue,Factor,Condition,Treatment,Replicate,bamReads,ControlID,bamControl,Peaks,PeakCallers
ws1,non,ATAC-seq,ws,,1,ATAC_seq_output/bwa/merged_library/ws_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP1.mLb.clN_peaks.broadPeak,bed
ws2,non,ATAC-seq,ws,,2,ATAC_seq_output/bwa/merged_library/ws_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP2.mLb.clN_peaks.broadPeak,bed
wf1,non,ATAC-seq,wf,,1,ATAC_seq_output/bwa/merged_library/wf_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/wf_REP1.mLb.clN_peaks.broadPeak,bed
wf2,non,ATAC-seq,wf,,2,ATAC_seq_output/bwa/merged_library/wf_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/wf_REP2.mLb.clN_peaks.broadPeak,bed
ms1,non,ATAC-seq,ms,,1,ATAC_seq_output/bwa/merged_library/ms_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ms_REP1.mLb.clN_peaks.broadPeak,bed
ms2,non,ATAC-seq,ms,,2,ATAC_seq_output/bwa/merged_library/ms_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ms_REP2.mLb.clN_peaks.broadPeak,bed
mf1,non,ATAC-seq,mf,,1,ATAC_seq_output/bwa/merged_library/ws_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP1.mLb.clN_peaks.broadPeak,bed
mf2,non,ATAC-seq,mf,,2,ATAC_seq_output/bwa/merged_library/ws_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP2.mLb.clN_peaks.broadPeak,bed

```

# Answer
> 2 votes
I have very little idea what's going on here, but in the absence of any other answers or help, I'll try to work through this from a no-knowledge start point to see if I can hopefully help shift this problem closer towards a solution.

The specific error you have encountered is `No valid contrasts/methods specified`. Looking at the DiffBind Vignette, it provides this information about defining contrasts:

> Before running the differential analysis, we need to tell DiffBind how to model the data, including which comparison(s) we are interested in. This is done using the `dba.contrast` function, as follows:

```
> tamoxifen <- dba.contrast(tamoxifen, reorderMeta=list(Condition="Responsive"))
> tamoxifen
11 Samples, 2845 sites in matrix:
       ID Tissue Factor  Condition  Treatment Replicate   Reads FRiP
 1 BT4741  BT474     ER  Resistant Full-Media         1  652697 0.16
 2 BT4742  BT474     ER  Resistant Full-Media         2  663370 0.15
 3  MCF71   MCF7     ER Responsive Full-Media         1  346429 0.31
 4  MCF72   MCF7     ER Responsive Full-Media         2  368052 0.19
 5  MCF73   MCF7     ER Responsive Full-Media         3  466273 0.25
 6  T47D1   T47D     ER Responsive Full-Media         1  399879 0.11
 7  T47D2   T47D     ER Responsive Full-Media         2 1475415 0.06
 8 MCF7r1   MCF7     ER  Resistant Full-Media         1  616630 0.22
 9 MCF7r2   MCF7     ER  Resistant Full-Media         2  593224 0.14
10  ZR751   ZR75     ER Responsive Full-Media         1  706836 0.33
11  ZR752   ZR75     ER Responsive Full-Media         2 2575408 0.22

```

After showing your contrast result, you explained a partial removal of the error:

> If I try to compare MS vs MF, I will get that error but if I remove MS and replace by other sample, I don't have that error anymore.

Specifically relating to MS vs MF, I can see why there would be an issue with the condition contrast, especially if it's using DESeq2 in the background. The short answer is that there is not enough variation within the groups.

DESeq2 cares about making sure that conditions added to the model are separable from samples. This is further explained in the DESeq2 manual here. The issue is that there's no way to distinguish between the control condition and the MS sample (or, similarly, between the diseased condition and the MF sample). In order to get past this error, you need both control and diseased in your MS group or in your MF group (or both, if you want to be *really* rigorous), and ideally a few more replicates of each condition.

---
Tags: r, bioconductor, atac-seq
---


title: Create bam from separate sam and header

text:
# Question
Title: Create bam from separate sam and header
If I have multiple SAM files without header, and a header saved in a separate file, how do I convert that to a BAM file?

One solution is to create intermediate files:

```
for samfile in *.sam
do
  cat my_header > tmp.sam
  cat $samfile >> tmp.sam
  samtools view -b tmp.sam > $samfile.bam
done

```

But that seems inelegant and I feel like I'm missing the better way.

# Answer
> 1 votes
I'd use GNU parallel. I've not tested this with real files but try it out, ensure the `{}` is passed to the process substitution:

```
mkdir with_header
parallel samtools view -b <(cat my_header {}) -o with_header/{} ::: *.bam

```

The command could be better - out file name could be done better but you get the concept.

---
Tags: samtools, bash
---


# Question
Title: Nextflow - Missing value declared as output parameter - output variable defined in script
I've normal-tumor human FASTQ file(s). Currently, I work with only one sample (both normal and tumor fastq R1 and R2), only.
I have alignment gatk mark duplicates, setupnmdtags, applybqsr, recalibrate steps before calling mutect2.

My input files are:
tiny\_tumor\_R1.fastq.gz
tiny\_tumor\_R2.fastq.gz
tiny\_normal\_R1.fastq.gz
tiny\_normal\_R2.fastq.gz

The expected output from the mutect2 is **tiny.mutect2out.vcf** that is ${sample\_name}".mutect2out.vcf"

I split `tiny_normal` using nextflow split and take first element of the split array, that is tiny.

However, I get error as:

> Caused by: Missing value declared as output parameter: sample\_name

I understand that the variable was defined in the script step. How do I retrieve at the output level?

Below is code of the mutect2\_process.nf

```
process mutect2 {

        maxForks 3
        debug true

 memory '15 GB'

        input:
        tuple val(sample_id_tumor),path(bqsrbam_tumor)
        tuple val(sample_id_normal),path(bqsrbam_normal)

        output:
        tuple val(sample_name),path(${sample_name}".mutect2out.vcf")

        script:
        def avail_mem = task.memory ? task.memory.toGiga() : 0
                def java_options = [
        avail_mem ? "-Xmx${avail_mem}G" : "",
        "-XX:+UseSerialGC",
    ]
        def sample_name=sample_id_normal.split('_normal')[0]

        """

        ml gatk/4.1.3.0
gatk --java-options "${java_options.join(' ')}" Mutect2 -L $params.target -I ${bqsrbam_tumor} -I ${bqsrbam_normal} -normal ${sample_name} -tumor ${sample_name} --output ${sample_name}".mutect2out.vcf" -R $params.hg38genome

        """
        }

```

# Answer
You would like to use the *sample\_name* variable outside the `script` block, but the `def` keyword defines a variable locally. All you'd need to do is remove the `def` keyword to make it work, but usually what you'd want to do is to define your *sample\_name* variable outside of your process block in your workflow block like in my previous answer. I would also call this variable something else like *pair\_id*, *donor\_name*, *case* or *patient* perhaps. I would say you actually have two sample names here: one tumor sample (i.e. *tiny\_tumor*) and one normal sample (i.e. *tiny\_normal*) from the same patient/donor/case.

You may also run into issues if *both* of your input channels are queue channels. You may not notice any issues if you're testing only a single pair of samples. Try just copying your tiny files a couple of times and giving your workflow a try, for example:

```
for i in tiny_{tumor,normal}_R{1,2}.fastq.gz ; do
    for j in foo bar baz; do
        cp "${i}" "${j}_${i#*_}"
    done
done

```

Most of the time what you'll want is one queue channel and one or more value channels when you require multiple input channels.
> 1 votes

---
Tags: nextflow
---


# Question
Title: Is there any function in Biopython to convert a DNA sequence from ambiguous to unambiguous?
I have a project to write in Python that requires to write a function that given degenerate DNA sequence (for example: KKGTACACCAG) sequence and a molecular weight interval, returns a list of all unambiguous sequences represented by the given sequence (encoding based on IUPAC alphabet: ***http://www.bioinformatics.org/sms/iupac.html***).

Is there any function in Biopython to convert a DNA sequence from degenerate to unambiguous base pairs?

Also, Is it possible to use the molecular weight function from Bio.SeqUtils library for an ambiguous DNA?

# Answer
> 3 votes
```
from Bio.Seq import Seq
coding_dna = Seq("AAYGANAGYCARAGYAAR")
print(coding_dna.translate())

```

This gives

> NXSQSK

So .translate() will read a degenerate code and provide the unambiguous amino acid, except where it is unable to do so e.g. GAN because this is a split between asp and glu. Thus GAR or GAY is fine, but GAN gives X

You would then use `itertools` by mapping (`map`) against the triplet code (`import Bio.Data.CodonTable`) to obtain all permutations. Doesn't work for X admittedly and stop codons a special case (3 possibilities).

`Bio.Data.CodonTable.list_possible_proteins` would solve the `X` problem, but stop codons remain a special case.

---

I seriously did like the Stackoverflow link in the comments, with the double map funtion ... very slick.

# Answer
> 3 votes
M\__'s approach is good, but, as explained, some subset of the unambiguous DNA sequences will not be recovered because of the ambiguity in reverse-translation due to the degeneracy of the genetic code. I don't know a `python` approach, but here's an `R` function that returns a vector of all the possible unambiguous DNA strings:

```
make_unambiguous <- function(dna) {
  require(tidyverse)
  require(S4Vectors)
  iupac <- tibble(code = c("A", "C", "G", "T",
                           "R", "Y", "S", "W",
                           "K", "M", "B", "D",
                           "H", "V", "N"),
                  base = c("A", "C", "G", "T",
                           "AG", "CT", "GC", "AT",
                           "GT", "AC", "CGT", "AGT",
                           "ACT", "ACG", "ACGT"))

  tibble(dna) |>
    separate_rows(dna, sep = '(?<=.)(?=.)') |>
    left_join(iupac, by = c("dna" = "code")) |>
    pull(base) |>
    str_split("") |>
    expand.grid(stringsAsFactors = FALSE) |>
    unite(col = sequence, sep = "") |>
    as_vector()
}

make_unambiguous(dna = "AAYGANAGYCARAGYAAR")

```

---
Tags: python, homework
---


# Question
Title: How do I delete an amino acid residue from a loop and then join the loose ends?
I want to delete an amino acid residue from a loop and then join the loose ends. I intend to increase the thermal stability of the protein by reducing the size of the loop and its agitation. Pymol just deletes the residue and breaks the loop.

I intend to increase the optimal temperature of the protein, even if it costs some reduction in activity. The objective is to reduce protein agitation due to temperature increase. I will try point mutations that increase hydrogen interactions and eventually delete a loop located glycine to reduce its length and consequently its thermal agitation. All this work is aimed at preserving the catalytic site, but some loss of activity is acceptable.

As a complement to my question, I present an excerpt from the paper by Mahdie Rahban et al (2022) ; Thermal stability enhancement: Fundamental concepts of protein engineering strategies to manipulate the flexible structure.

> Loop deletion or shortening by insertion of proline or deletion of Gly, coupled with loop anchoring by increasing the number of salt bridges, hydrogen bonds, or hydrophobicity, are good strategies to enhance the overall rigidity of a protein conformation leading to better compactness and a decreased number of internal cavities.

# Answer
> 0 votes
It depends on what software you have in hand and want to use. But basically, when you delete/add a residue, you can't just stick two parts together since it affects the whole conformation.

I suggest using AlphaFold2, and input it with different sequences (including deletion mutation) to see how it affects the structure.

Here is a link to AF2 colab: https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb

---
Tags: proteins, protein-structure
---


# Question
Title: Understanding DESeq2 output for resultNames
I'm following along with some pre written code to learn how to do differential expression. I am using the same data set and have gone back through the code and copy and pasted what they used and the follow along gets the two headings rather than the single one that I am getting, so I'm a bit confused how I'm getting different results.

I have some data which I have put through DESeq but the output for resultsNames isn't what I expect. The code is

```
dds <- DESeqDataSetFromMatrix(countData = exprs(bottomly.eset), colData = pData(bottomly.eset), design = ~ strain)
dds <- DESeq(dds)

```

but this code produces this output `"Intercept" "strain_DBA.2J_vs_C57BL.6J"` whereas the output that I intended was `"Intercept" "strain_DBA.2J" "C57BL.6J"`

Is there a simple way to change this within the code that I have or is this something related to the dataset I used?

# Answer
> 2 votes
DESeq2 is a package designed for calculating *differential* expression, in other words for comparing expression between two states. The primary results it produces are therefore comparative values (e.g. comparing DBA.2J and C57BL.6J in your case).

Given that this does not seem to be your intended calculation, it would be helpful if you could provide more context around what you want to use DESeq2 for. If it's not for differential expression, there's likely a better tool / approach for solving your problem.

# Answer
> 2 votes
So I was following along with a coding lesson from Youtube learning how to do differential expression in R, if the data had been my own I probably would have been able to figure out what needed to be done. The code that worked for me was

`dds <- DESeqDataSetFromMatrix(countData = exprs(bottomly.eset), colData = pData(bottomly.eset), design = ~ strain)`

`dds <- DESeq(dds, betaPrior = TRUE)`

I was then able to use this code below which was resulting in the error previously `contrast.deseq2 <- list(c("strainC57BL.6J", "strainDBA.2J"))`

Thank you to everyone who commented on this question

---
Tags: differential-expression, deseq2
---


# Question
Title: Sql query slow on Chembl postgres database
Doing complex sql query on Chembl database and it takes 40 minutes on the machine from lab. I was thinking using dask, extract each table and export to parquet, when do the data processing . Is this a good approach or there is a better alternative?

# Answer
> 0 votes
I suspect the issue is you want to parallelise (multi-threading really) the job to reduce run time.

`parquet` is recommended. It is simply an archive format that allows compression and is often used with a `pandas` dataframe. Use `parquet` to store data by all means, but thats easy via `pyarrow` and the `pandas` command

```
    df.to_parquet('/User/OP/location/file,parquet' ,engine='pyarrow', compression='gzip')

```

You will need to install `pyarrow` via pip or conda, there is a faster version (forget its name, something like `fastparquet` or something used by default). Yeah it's a very good archiving format, in fact its industrial standard and I should use it more.

Python has loads of ways to parallelise of which dask is one way. dask is an entire framework. To leverage dask to parallelise postgres I would refer to this post.

I don't use dask, I am aware it is a solution in high performance Python. Personally I would use this approach directly via `multiprocessor`

```
from psycopg2 import connect
import multiprocessor

```

... etc ...

I would suggest its horses for courses, I simply know how to multiprocess in Python so I don't bother with dask.

To answer the question in summary

* parquet is the best, if SPEED is not important. If speed is important then `feather`
* dask is mainstream within high performance Python and its doable. There's no reason not to
* There are other multithread Postgress approaches and thats what I would do, but if you just getting into parallelisation you may as well use dask. It will take some learning though.

---

If you're not into parallelisation and you're not writing a data pipeline for publication I would simply consider getting a machine with more RAM. It shouldn't be taking that length of time, parallelisation will solve the run time - no question - but it's the time it will take to get to grips with a parallel solution. What is likely occurring is the job and database size is exceeding the RAM capacity of the machine. At a wild random guess a 16 or 32 Meg RAM chip would collapse the runtime to a few minutes and thats what you are looking for.

---

From the comments.

`parquet` is really for data compression and archiving. The issue is I/o is reasonably fast but not lightening fast. It is probably slower than `pickle` (it will be slower if you're compressing via gzip as well).

For a 'temporary' save (couple months) and really fast i/o then its `feather`. There is no dispute about that.

If you are keeping your data a long time and need to save space its `parquet`. Performance-wise `feather` is way faster than `parquet`. I've never used `feather` BTW I use `pickle` or `parquet`

To give some idea a 500Meg flat file compressed an archived via `parquet` resulting in a 50Meg file, and *will take several minutes to write/import*. Thats 'cause it's doing loads of stuff. My suspicion is that this too slow for you.

---
Tags: python, database, data-retrieval, cheminformatics
---


# Question
Title: How to open .weir.fst files witout losing data?
I have .vcf file with over 78 000 variations in 95 experiment subjects. I used vcftools to extract fixation index score for general differentiation between Case and Control groups.

```
./vcftools --vcf SCH_freebayes_sort_ostE_SnpEff_SnpSift.vcf --weir-fst-pop Case.txt --weir-fst-pop Control.txt --out FstExon_Case_vs_Control

```

```
VCFtools - 0.1.17
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
    --vcf SCH_freebayes_sort_ostE_SnpEff_SnpSift.vcf
    --weir-fst-pop Case.txt
    --weir-fst-pop Control.txt
    --keep Case.txt
    --keep Control.txt
    --out Fst_TestData_nochange

Warning: Expected at least 2 parts in INFO entry: ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
Warning: Expected at least 2 parts in INFO entry: ID=PRO,Number=1,Type=Float,Description="Reference allele observation count, with partial observations recorded fractionally">
Warning: Expected at least 2 parts in INFO entry: ID=PAO,Number=A,Type=Float,Description="Alternate allele observations, with partial observations recorded fractionally">
Warning: Expected at least 2 parts in INFO entry: ID=SRP,Number=1,Type=Float,Description="Strand balance probability for the reference allele: Phred-scaled upper-bounds estimate of the probability of observing the deviation between SRF and SRR given E(SRF/SRR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=SAP,Number=A,Type=Float,Description="Strand balance probability for the alternate allele: Phred-scaled upper-bounds estimate of the probability of observing the deviation between SAF and SAR given E(SAF/SAR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=AB,Number=A,Type=Float,Description="Allele balance at heterozygous sites: a number between 0 and 1 representing the ratio of reads showing the reference allele to all reads, considering only reads from individuals called as heterozygous">
Warning: Expected at least 2 parts in INFO entry: ID=ABP,Number=A,Type=Float,Description="Allele balance probability at heterozygous sites: Phred-scaled upper-bounds estimate of the probability of observing the deviation between ABR and ABA given E(ABR/ABA) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=RPP,Number=A,Type=Float,Description="Read Placement Probability: Phred-scaled upper-bounds estimate of the probability of observing the deviation between RPL and RPR given E(RPL/RPR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=RPPR,Number=1,Type=Float,Description="Read Placement Probability for reference observations: Phred-scaled upper-bounds estimate of the probability of observing the deviation between RPL and RPR given E(RPL/RPR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=EPP,Number=A,Type=Float,Description="End Placement Probability: Phred-scaled upper-bounds estimate of the probability of observing the deviation between EL and ER given E(EL/ER) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=EPPR,Number=1,Type=Float,Description="End Placement Probability for reference observations: Phred-scaled upper-bounds estimate of the probability of observing the deviation between EL and ER given E(EL/ER) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=CIGAR,Number=A,Type=String,Description="The extended CIGAR representation of each alternate allele, with the exception that '=' is replaced by 'M' to ease VCF parsing.  Note that INDEL alleles do not have the first matched base (which is provided by default, per the spec) referred to by the CIGAR.">
Warning: Expected at least 2 parts in FORMAT entry: ID=GQ,Number=1,Type=Float,Description="Genotype Quality, the Phred-scaled marginal (or unconditional) probability of the called genotype">
Warning: Expected at least 2 parts in FORMAT entry: ID=GL,Number=G,Type=Float,Description="Genotype Likelihood, log10-scaled likelihoods of the data given the called genotype for each possible genotype generated from the reference and alternate alleles given the sample ploidy">
Warning: Expected at least 2 parts in INFO entry: ID=WTD,Number=0,Type=Flag,Description="Is Withdrawn by submitter If one member ss is withdrawn by submitter, then this bit is set.  If all member ss' are withdrawn, then the rs is deleted to SNPHistory">
Warning: Expected at least 2 parts in INFO entry: ID=WTD,Number=0,Type=Flag,Description="Is Withdrawn by submitter If one member ss is withdrawn by submitter, then this bit is set.  If all member ss' are withdrawn, then the rs is deleted to SNPHistory">
Warning: Expected at least 2 parts in INFO entry: ID=TOPMED,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on TOPMed, starting with the reference allele followed by alternate alleles as ordered in the ALT column. The TOPMed minor allele is the second largest value in the list.">
Warning: Expected at least 2 parts in INFO entry: ID=TOPMED,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on TOPMed, starting with the reference allele followed by alternate alleles as ordered in the ALT column. The TOPMed minor allele is the second largest value in the list.">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=MTP,Number=0,Type=Flag,Description="Microattribution/third-party annotation(TPA:GWAS,PAGE)">
Warning: Expected at least 2 parts in INFO entry: ID=NOC,Number=0,Type=Flag,Description="Contig allele not present in variant allele list. The reference sequence allele at the mapped position is not present in the variant allele list, adjusted for orientation.">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=MUT,Number=0,Type=Flag,Description="Is mutation (journal citation, explicit fact): a low frequency variation that is cited in journal and other reputable sources">
Warning: Expected at least 2 parts in INFO entry: ID=SAO,Number=1,Type=Integer,Description="Variant Allele Origin: 0 - unspecified, 1 - Germline, 2 - Somatic, 3 - Both">
Warning: Expected at least 2 parts in INFO entry: ID=SAO,Number=1,Type=Integer,Description="Variant Allele Origin: 0 - unspecified, 1 - Germline, 2 - Somatic, 3 - Both">
Warning: Expected at least 2 parts in INFO entry: ID=SAO,Number=1,Type=Integer,Description="Variant Allele Origin: 0 - unspecified, 1 - Germline, 2 - Somatic, 3 - Both">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=PM,Number=0,Type=Flag,Description="Variant is Precious(Clinical,Pubmed Cited)">
Warning: Expected at least 2 parts in INFO entry: ID=AF_popmax,Number=A,Type=Float,Description="Maximum allele frequency across populations (excluding samples of Ashkenazi, Finnish, and indeterminate ancestry)">
Warning: Expected at least 2 parts in INFO entry: ID=AF_popmax,Number=A,Type=Float,Description="Maximum allele frequency across populations (excluding samples of Ashkenazi, Finnish, and indeterminate ancestry)">
Keeping individuals in 'keep' list
After filtering, kept 95 out of 95 Individuals
Outputting Weir and Cockerham Fst estimates.
    Fst: Only using diploid sites.
Weir and Cockerham mean Fst estimate: -0.00031927
Weir and Cockerham weighted Fst estimate: -0.00047509
After filtering, kept 78758 out of a possible 78758 Sites
Run Time = 4.00 seconds

```

But when I tried to open .weir.fst files with notepad or excel, I lost over 9 000 rows of data. I am sure that there is problem with how I am opening the .weir.fst files.

Please, how are these files opened properly?

# Answer
> 0 votes
Oh this is fairly easy. The file size has exceed the capacity of Excel (50000 rows??). This doesn't affect downstream analysis because the rows are there you just can't visualise them.

**BBEdit** will open it if this is Mac (free from AppStore - really good BTW). If this is Linux then `vi` or `vim`. A Windows version of BBEdit is **Notepad++**, please not Notepad.

You can request the number of lines (bash Linux) via

```
wc weir.fst

```

This will return the number of lines in the file, if it's \>50000, i.e. beyond the number of the Excel limit, that's your answer. You just need to find an alternative viewer

---
Tags: vcftools, populations, population-genetics, fst
---


# Question
Title: Choosing the best number of species for assay
I am doing amplicon sequencing of a virus across many different regions. Lets say I have 20k unique variants of the same virus that I put into my pcr assay and after sequencing and amplifications I am left with, 19k variants that appear. But many of the variants appear at really low read count and I don't know if that read is real or just noise. I have a threshold T that I say if a variant appears above T I'm counting it as real. Out of the 20k variants about only 15% of variants appear to be "real". I have some data on different input species counts (2k, 10k, 20k, 40k), and their corresponding fraction of "real".

My question is how do I determine which input variants count is best? I only looked at 4 different values and maybe the correct value isn't sampled, I obviously can't try everything between 2k and 40k. I think I want a balance between largest fraction of "real" variants and total number of variants.

Is there any research (or better yet code) that answers/discusses this problem?

Thanks in advanced.

# Answer
> 1 votes
The underlying answer to the question is have you experimentally controlled for *in vitro* artefacts? Thats taking a cloned virus and performing the same RT-PCR. If not thats makes things very difficult. If there isn't a cloned virus thats a problem. You might need to ask your wet-lab colleagues what old school "cloning a virus" means, or Biology SE will explain this. The reason this is important is given below.

The first step to the question is simply to plot a histogram against copy number. What you will find is a heavily skewed distribution where the copy number of one variant will represent \>40% (at a wild guess) of all the viral copy number. Moreover, a handful (less than 100) of variants will comprise 90% of all the copy number. Once that is established it is much easier to assess what is a true variant. What I'm saying is the difference is one log minimum between what you think you are looking at and what is there.

The sources of artefact:

* RT-PCR artefacts, particularly if there is an RT step: Taq polymerase alone has an error rate of 1 in 1000 and RT has a very high error rate, worse it's the first step.

What DNA polymerase was used?

The key issue is separating the RT-PCR artefacts from naturally occurring RNA polymerase errors *in situ* is not possible without a control. If a frequency analysis cut-off strategy is used it needs to be an ultra-conservative threshold, i.e. I mean ridiculously high and this means lot of variants are going to be dumped.

**Why?** Again it's back to the problem whether the variants were generated *in vitro*, i.e. artefacts, or *in situ*, i.e. resulting from RNA polymerase. Generating a cut-off within a skewed distribution is not trivial *per se* but extremely difficult when the ratio between artefact and genuine variant has not been controlled.

---

There appears some sort of control. Maximising the threshold: thats a function of the control, so the control needs to accommodate a model of the different numbers of input variables; it will not be a linear function. Basically what input will give me the best representation of non-artifactual variants, that is what the control has to find out. Thus, for every input, you want a control.

Again, if the control is lacking, in any context, then the rule of thumb is to be ultra-conservative, not trying to maximise the number of variants, but minimise them, particularly if the control isn't perfect. Most virological reviewers raise concern with the numbers been reported here on the information given (obviously we don't know the virus).

* NGS sequencing is another source of artefact - I should have mentioned there are some famous disasters here. I'm sure you NGS is really good, but a small number of investigators have deposited a disproportionate number of artefact ridden sequences - that everyone knows about.

---

from comments (again)

A control of 4 is tiny in data science especially a GLM approach. However, other studies would have been performed - possibly lots of them if this was e.g. HIV - and controls would have been run. There are two major strategies (really 3) for biologically cloning a virus. Really you'll need 100+ controls, that may not be feasible but as many within 4 to 100 as are available. GLM would be a good approach IMO. It would need stating that this does not truly control the rate of artefacts, but it's much closer.

---

Final comments ... don't want to assess thresholds.

I gave this question some thought and I think it would work via machine learning (what I specialise in). However, it may be difficult to find appropriate published controls and/or augmentation. Personally, I'm satisfied this would work if done correctly, but implementation is fairly advance ML. Given the amount of time I've spent on thinking this question through I wouldn't use an alternative approach, because I think it would be valid. Augmentation BTW is a generic term in machine learning and deep learning to increase the size of the training set.

---
Tags: ngs, variant-calling, virus, pcr, qpcr
---


# Question
Title: selenium webdriver is stuck after downloading and saving a pdf file from pubmed central website
I'm trying to download pdf from pubmed central using selenium webdriver in python. The code is provided below which downloads the pdf fine but, it does not quit automatically. i.e. it is sort of stuck executing `browser.get(href)` line in `download_pmc_pdf()`. What could be the reason?

```
import os
from selenium import webdriver
from selenium.webdriver.common.by import By
from selenium.webdriver.firefox.options import Options
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC
ncbi_url = "https://www.ncbi.nlm.nih.gov/"

download_dir = "pdf_downloaded"

if not os.path.exists(download_dir):
    os.mkdir(download_dir)

options = Options()
options.set_preference('permissions.default.stylesheet', 2)
options.set_preference('permissions.default.image', 2)
options.set_preference('browser.download.folderList', 2)
options.set_preference('browser.download.dir', os.path.abspath(download_dir))
options.set_preference('browser.download.manager.showWhenStarting', False)
options.set_preference('browser.helperApps.neverAsk.saveToDisk', 'application/pdf')
options.set_preference('pdfjs.disabled', True)

browser = webdriver.Firefox(options=options)

###############################################################################
def download_pmc_pdf(pmcid_key):
    pmc_url = ncbi_url + "pmc/articles/" + pmcid_key + "/"
    browser.get(pmc_url)
    try:
        wait = WebDriverWait(browser, 10)
        pdf_links = wait.until(EC.presence_of_all_elements_located((By.CSS_SELECTOR, '.pmc-sidebar__formats li a')))

        # Iterate over all links and navigate to the link containing '.pdf' in the href attribute
        for link in pdf_links:
            href = link.get_attribute('href')
            if '.pdf' in href:
                browser.get(href)
                break  # exit the loop as we have found the PDF link
        print("line beofre browser.quit()")
        browser.quit()  # quit the browser outside the loop
    except Exception as e:
        print(e)

###############################################################################
pmcid = "2788555"
download_pmc_pdf(pmcid)
browser.quit()

```

# Answer
> 2 votes
I had to execute a JavaScript code that opens the pdf link in a new window/tab. After that, just closed this newly opened window/tab and continued with the script execution. Here is the updated code:

```
import os
import time
import csv

from selenium import webdriver
from selenium.webdriver.common.by import By
from selenium.webdriver.firefox.options import Options
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC

ncbi_url = "https://www.ncbi.nlm.nih.gov/"

download_dir = "pdf_downloaded"

if not os.path.exists(download_dir):
    os.mkdir(download_dir)

options = Options()
options.set_preference('permissions.default.stylesheet', 2)
options.set_preference('permissions.default.image', 2)
options.set_preference('browser.download.folderList', 2)
options.set_preference('browser.download.dir', os.path.abspath(download_dir))
options.set_preference('browser.download.manager.showWhenStarting', False)
options.set_preference('browser.helperApps.neverAsk.saveToDisk', 'application/pdf')
options.set_preference('pdfjs.disabled', True)

###############################################################################
def download_pmc_pdf(pmcid_key):
    browser = webdriver.Firefox(options=options)
    pmc_url = ncbi_url + "pmc/articles/" + pmcid_key + "/"
    browser.get(pmc_url)
    try:
        wait = WebDriverWait(browser, 10)
        pdf_links = wait.until(EC.presence_of_all_elements_located((By.CSS_SELECTOR, '.pmc-sidebar__formats li a')))

        # Iterate over all links and navigate to the link containing '.pdf' in the href attribute
        for link in pdf_links:
            href = link.get_attribute('href')
            if '.pdf' in href:
                # Open the pdf link in a new window
                browser.execute_script(f'window.open("{href}","_blank");')
                break

        time.sleep(15)  # delay for 10 seconds
        # After the file is downloaded:
        with open(f'{download_dir}/file_mapping.csv', 'a', newline='') as file:
            writer = csv.writer(file)
            writer.writerow([pmcid_key, os.path.basename(href)])
        ## Switch to the new window and close it
        ## (this is a hack to deal with bowser freezing upon downloading)
        browser.switch_to.window(browser.window_handles[-1])
        browser.close()

    except Exception as e:
        print(e)
    finally:
        browser.quit()  # quit the browser outside the loop
###############################################################################
pmcid = "2788555"
download_pmc_pdf(pmcid)

```

# Answer
> 1 votes
Pubmed Central stores pdfs for all available open access articles on a public-facing ftp server (https://ftp.ncbi.nlm.nih.gov/pub/pmc/oa\_pdf/). You can query the database using the OA Web Service API, e.g. with cURL

```
# for pubmed id "10026098"
curl https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=10026098
<OA><responseDate>2023-07-17 01:39:13</responseDate><request id="PMC10026098">https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=10026098</request><records returned-count="2" total-count="2"><record id="PMC10026098" citation="Ther Adv Med Oncol. 2023 Mar 18; 15:17588359231161409" license="CC BY-NC" retracted="no"><link format="tgz" updated="2023-03-23 12:10:30" href="ftp://ftp.ncbi.nlm.nih.gov/pub/pmc/oa_package/83/eb/PMC10026098.tar.gz" /><link format="pdf" updated="2023-03-20 09:34:33" href="ftp://ftp.ncbi.nlm.nih.gov/pub/pmc/oa_pdf/83/eb/10.1177_17588359231161409.PMC10026098.pdf" /></record></records></OA>

# get the pdf link from the API (ftp storage)
curl https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=10026098 | sed -n 's/.*\"\(ftp.*.pdf\)\".*/\1/p'
ftp://ftp.ncbi.nlm.nih.gov/pub/pmc/oa_pdf/83/eb/10.1177_17588359231161409.PMC10026098.pdf

# download the pdf
pdf="$(curl https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=10026098 | sed -n 's/.*\"\(ftp.*.pdf\)\".*/\1/p')"
curl -O "$f"

# download many pdfs from a list of ids
cat list_of_ids.txt
10026098
9341183

while read -r id_from_list
do
    pdf="$(curl https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=$id_from_list | sed -n 's/.*\"\(ftp.*.pdf\)\".*/\1/p')"
    curl -O "$pdf"
done < list_of_ids.txt

```

This approach is preferable to selenium scraping for many reasons, and if you have any problems the NCBI provides documentation (https://www.ncbi.nlm.nih.gov/pmc/tools/ftp/#indart) and offers support (you can email them at pubmedcentral@ncbi.nlm.nih.gov).

---
Tags: python
---


# Question
Title: how to merge more than two sample in Seurat?
I would like to merge more than two sample in the Seurat, and the mergeseurat can only merge two sample. So what should I do now. The screenshot is my script.

# Answer
> 4 votes
You have fed arguments to the `MergeSeurat()` function that it does not expect. In terms of objects, `MergeSeurat()` accepts only 2 arguments `object1` and `object2`, please do `?MergeSeurat` and see yourself.

Seurat is probably the best documented single cell package, your exact question has its own dedicated vignette: https://satijalab.org/seurat/v3.0/merge\_vignette.html

As the vignette above is for Seurat v3 and in case you would like to use Seurat v2: You can create a named list of Seurat objects to be merged and then use `reduce()` of the `purrr` package.

```
# create a named list of seurat objects to be merged
# code not shown

# optional but probably a good idea
# rename cells using object names as prefix
for (i in names(seurat_object_list)) {
  seurat_object_list[[i]] <- RenameCells(seurat_object_list[[i]],
                                         add.cell.id = i)
}

# merge all the objects in the list
merged_combined <- reduce(seurat_object_list,
                          MergeSeurat,
                          do.normalize = FALSE)

```

# Answer
> 2 votes
The `MergeSeurat` command is from Seurat v2.4 and only accepts two objects as parameters. If you need to merge more than one you can first merge two, then merge the combined object with the third and so on.

In the newer Seurat v3.0 this is replaced by the `merge` command that can have a named list of Seurat objects as input

```
# merge two objects
merge(x = pbmc_small, y = pbmc_small)
# to merge more than two objects, pass one to x and a list of objects to y
merge(x = pbmc_small, y = c(pbmc_small, pbmc_small))

```

If you want to perform batch correction at the same time as merging the data you should have a look at the `IntegrateData` functionality of Seurat v3.0, explained in the corresponding vignettes:

https://satijalab.org/seurat/v3.0/immune\_alignment.html
https://satijalab.org/seurat/v3.0/integration.html

A similar workflow for Seurat v2.4 can be cound here:

https://satijalab.org/seurat/v2.4/immune\_alignment.html

# Answer
> 2 votes
As of 2023 (or Seurat \>4 or \>3?), simply use `merge` and write

```
combined.obj <- merge(x = ls.Seurat[[1]], y = ls.Seurat[2:l(ls.Seurat)], add.cell.ids = vec_sample_IDs)

```

Since `y` can take a list.

# Answer
> 1 votes
For now, Seurat does not have an option to merge more than two samples at a time. You can write a function to run the `MergeSeurat` function in a loop. Following issue in Github has an example function for multiple samples.

https://github.com/satijalab/seurat/issues/432

---
Tags: rna-seq, seurat, single-cell
---


title: nextflow: unable to find index

text:
# Question
Title: nextflow: unable to find index
I've a pipeline with many processes.
One of the steps is to run lancet tool. Lancet tool runs per chromosome (1-22, x,y). The input in the lancet step is output from the applybqsr (gatk).

There are tumor and normal bam generated along with the index. human genome FASTA has a index. When I run lancet through the main.nf I run into an error

> ERROR: index not found for BAM file

I've a config file as:

```
params {

        hg38genome="/sc/arion/reference_genome/human/hg38/Homo_sapiens_assembly38.fasta.gz"
        outdir = '/sc/arion/output_pipeline_nextflow/test_tiny_datasets/'
        dbsnp="/sc/arion/projects/Daphni2_Pipeline_resources/resourceDBsnp146/dbsnp_146.hg38.vcf.gz"
        target="/sc/arion/projects/Twist_Exome_Target_hg38.padded.interval_list.bed/Twist_Exome_Target_hg38.padded.interval_list.bed"

}
        withName: applybqsr {
        publishDir = [
                [
                path: "${params.outdir}/applybqsr/normal/",
                mode: 'copy',
                pattern:'*normal*'
                ],
                [
                path: "${params.outdir}/applybqsr/tumor/",
                mode: 'copy',
                pattern:'*tumor*'
                ]
        ]
        }
        //applybqsr ends

        withName: lancet {
        publishDir = [
                path: "${params.outdir}/lancet/",
                mode: 'copy'
        ]
        }
        //lancet ends

```

I have applybqsr as:

```
process applybqsr{

        maxForks 3
        debug true

        memory '6 GB'

        input:
        tuple val(sample_id_tumor), path(markbam_tumor)
        tuple val(sample_id_tumor), path(recaltable_tumor)
        tuple val(sample_id_normal), path(markbam_normal)
        tuple val(sample_id_normal), path(recaltable_normal)

        output:
        tuple val(sample_id_tumor), path("${sample_id_tumor}.sorted.markdup.recal.bam"), emit: recal_bam_tumor
        tuple val(sample_id_tumor), path("${sample_id_tumor}.sorted.markdup.recal.bai"),  emit: recal_bai_tumor
        tuple val(sample_id_normal), path("${sample_id_normal}.sorted.markdup.recal.bam"), emit: recal_bam_normal
        tuple val(sample_id_normal), path("${sample_id_normal}.sorted.markdup.recal.bai"),  emit: recal_bai_normal

        script:
        """
        ml gatk/4.1.3.0 samtools/1.9
        samtools index ${markbam_tumor}
        samtools index ${markbam_normal}

        gatk ApplyBQSR \\
        -R $params.hg38genome \\
        -I ${markbam_tumor} \\
        --bqsr-recal-file ${recaltable_tumor} \\
        -O ${sample_id_tumor}.sorted.markdup.recal.bam \\
        -L $params.target \\
        --create-output-bam-index true \\
        --java-options -Xmx${task.memory.toGiga()}g

        gatk ApplyBQSR \\
        -R $params.hg38genome \\
        -I ${markbam_normal} \\
        --bqsr-recal-file ${recaltable_normal} \\
        -O ${sample_id_normal}.sorted.markdup.recal.bam \\
        -L $params.target \\
        --create-output-bam-index true \\
        --java-options -Xmx${task.memory.toGiga()}g

        """

}
workflow.onComplete {
        log.info ( workflow.success ? "Done applyBQSR!" : "Oops .. something went wrong in applyBQSR" )
}

```

I've lancet.nf as:

```
process lancet {
        maxForks 3
        debug true

        input:
        tuple val(sample_id_tumor),path(bqsrbam_tumor)
        tuple val(sample_id_normal),path(bqsrbam_normal)
        val chromosome

        output:
        path "CHR${chromosome}_${sample_name}.lancet.vcf"

        script:

        sample_name=sample_id_normal.split('_normal')[0]

"""
ml lancet

lancet --tumor ${bqsrbam_tumor} --normal ${bqsrbam_normal} --ref $params.hg38genome --reg "chr"$chromosome --num-threads 8 > "CHR"${chromosome}"_${sample_name}.lancet.vcf"
"""

}

workflow.onComplete {
          log.info(workflow.success ? "Done lancet!" : "Oops .. something went wrong in lancet" )
 }

```

I have main.nf as:

```
params.outdir_index="/sc/arion/output_pipeline_nextflow/test_tiny_datasets/bwa_index"
params.raw_files_normal="/sc/arion/tiny/all_tumors_normals/*_normal_*{R1,R2}.fastq"
params.raw_files_tumor="/sc/arion/tiny/all_tumors_normals/*_tumor_*{R1,R2}.fastq"

include { FASTP} from './fastp_process.nf'
include {bwa_index} from './index_process.nf'

include { align_bwa_mem} from './bwamem_process_already_index.nf'
include { gatk_markduplicates} from './gatk_markduplicates_process.nf'
include {setupnmdtags} from './setupnmdtags_process.nf'
include { recalibrate_bam } from './recalibratebam_process.nf'
include { applybqsr } from './applybqsr_process.nf'
include { mutect2 } from './mutect2_process.nf'
include { lancet } from './lancet_process.nf'

workflow {

        normal_read_pairs=Channel.fromFilePairs(params.raw_files_normal)
        tumor_read_pairs=Channel.fromFilePairs(params.raw_files_tumor)
        //FASTP   align_bwa_mem        gatk_markduplicates         setupnmdtags        recalibrate_bam

        applybqsr(gatk_markduplicates.out.sorted_bam_mark_duplicates_tumor,recalibrate_bam.out.recalib_table_tumor,gatk_markduplicates.out.sorted_bam_mark_duplicates_normal,recalibrate_bam.out.recalib_table_normal)

        mutect2(applybqsr.out.recal_bam_tumor,applybqsr.out.recal_bam_normal)

        chr=Channel.of( 1..22, 'X', 'Y' )
        lancet(applybqsr.out.recal_bam_tumor,applybqsr.out.recal_bam_normal,chr)
}

```

Lancet runs fine if I run it hard coded into a process code block with same input (tumor-normal bam) and human genome FASTA. How do I enable index files to be staged into the lancet process generated one step before?

# Answer
I am reluctant to go through this wall of code, but generally it is advisable to always keep bam files and its index in a tuple so the index is staged with its bam file.

For example a process that creates the bam and now emits a tuple with bam and bai file:

```
(...)
output:
tuple path("${id}.bam"), path("${id}.bam.bai"), emit: tuple_bam

```

and then a process that receives such a `tuple_bam`:

```
(...)
input:
tuple path(bam), path(bai)

```

The first entry is the bam, the second the index.
> 2 votes

---
Tags: nextflow
---


# Question
Title: Assigning multiple GO IDs to an ID: Should I split or use the first or last ID in the hierarchy?
I have a dataset that contains IDs and associated Gene Ontology (GO) IDs. Some of the IDs have multiple GO IDs assigned to them, and I'm unsure how to handle this situation. I would like to determine whether Should I split all Go id or I should use the first or last GO ID in the hierarchy when assigning the GO ID to the ID.

Here is an example of my data:

```
3667 FEHLGCFF_05007                                             GO:0006637|GO:0047617
3673 FEHLGCFF_05919                                  GO:0016020|GO:0022857|GO:0055085
3684 FEHLGCFF_00369                                             GO:0022857|GO:0055085
3695 FEHLGCFF_00369                                  GO:0016020|GO:0022857|GO:0055085
3696 FEHLGCFF_00369                                                        GO:0022857
3702 FEHLGCFF_06641                                                        GO:0016627
3707 FEHLGCFF_06641                                                        GO:0003995
3721 FEHLGCFF_04414                                                        GO:0046872
3725 FEHLGCFF_02968                                  GO:0003700|GO:0006352|GO:0006355

```

How should I handle cases where an ID has multiple GO IDs assigned to it? Should I split all Go id or Should I use the first GO ID, the last GO ID, or is there another approach I should consider?

Any guidance or best practices would be greatly appreciated. Thank you!

Here is My code using R, it it Ok

```
library(tidyverse)

# Read the go-basic.obo file
go_data <- readLines("go-basic.obo")

# Extract the GO ID and Gene ID from the unique_data dataframe
unique_data <- unique_data %>%
  separate_rows(GO_ID, sep = "\\|") %>%
  select(V1, GO_ID) # Replace with the correct column names from your unique_data dataframe

# Create an empty dataframe to store the results
result <- data.frame(Gene_ID = character(), GO_ID = character(), Name = character(), Namespace = character(), stringsAsFactors = FALSE)

# Process each line in the go-basic.obo file
for (i in 1:(length(go_data) - 1)) {
  line <- go_data[i]

  if (grepl("^id:", line)) {
    go_id <- strsplit(line, " ")[[1]][2]

    # Check if the current GO ID exists in the unique_data dataframe
    if (go_id %in% unique_data$GO_ID) {
      # Filter the unique_data dataframe for the current GO ID
      filtered_data <- unique_data %>% filter(GO_ID == go_id)

      # Extract the name and namespace information
      name <- ""
      namespace <- ""

      while (!grepl("^is_a:", go_data[i])) {
        if (grepl("^name:", go_data[i])) {
          name <- strsplit(go_data[i], "name: ")[[1]][2]
        } else if (grepl("^namespace:", go_data[i])) {
          namespace <- strsplit(go_data[i], "namespace: ")[[1]][2]
        }
        i <- i + 1
      }

      # Add the information to the result dataframe
      temp_df <- data.frame(Gene_ID = filtered_data$V1, GO_ID = go_id, Name = name, Namespace = namespace, stringsAsFactors = FALSE) # Replace with the correct column names
      result <- bind_rows(result, temp_df)
    }
  }
}

# Print the result dataframe
result
nrow(result)

```

# Answer
> 1 votes
You might want to look at Revigo to simplify multiple GO terms into their parent terms. I've used rrvgo extensively to simplify enrichment results - this is the R implementation of Revigo.

It's not clear how exactly you are obtaining these GO assignments. If these are the significant hits for each feature, I'd take the union of all significant GO terms (i.e. everything in your second column combined) and pass it to rrvgo. Then, you can select the parent term associated with each row. I don't guarantee each row will be collapsed to a single term, but it should make whatever you are doing downstream much more tractable.

---
Tags: r, gene-ontology
---


# Question
Title: What kind of data do I need to be able to complete Differential Expression using Limma
I am very new to bioinformatics but have some background in coding with R. I have been asked to find a data set and complete differential expression using the Limma package in R. I have tried a few data sets but haven't been able to get any code to work.

What kind of data do I need to have in order for Differential Expression to work with Limma R? I have tried load and read.table, is there a read function that I am missing?

I was advised to look on Gene Expression Omnibus to find published data sets to work on but again I am struggling to find any that will work.

Thank you

# Answer
> 0 votes
Limma is relatively generic in what it accepts. Typically it is either data on log2 scale that are either already normalized or need some of the inbuilt normalization methods (such as quantile normalization), or in case of RNA-seq you can provide raw counts which then can run through the limma-voom pipeline.

The format is always a matrix of counts/values, with genes/observations in rows and samples/specimen in columns.

For practice, you could follow the RNA-seq/limma-voom section (chapter 15) of the limma manual which has a RNA-seq case study available at https://subread.sourceforge.net/RNAseqCaseStudy.html with example data and code.

GEO has a lot of published data, often microarray gene expression values or RNA-seq, but lets discuss this in a separate question if you need it to keep this here focused on limma.

---
Tags: r, differential-expression
---


# Question
Title: nextflow: merge lane files from samples forward and reverse read
I'd like to merge files from different lanes for forward and reverse reads. For e.g., if a sample A has three lanes in forward, and three reverse, ultimately one reverse and one forward is to be generated. I have normal and tumor for each sample. N as normal and T as tumor.
The lanes files are for e.g. as:

```
MM-265-DNA-N-01-01_L003_R1_001.fastq.gz
MM-265-DNA-N-01-01_L003_R2_001.fastq.gz
MM-265-DNA-N-01-01_L004_R1_001.fastq.gz
MM-265-DNA-N-01-01_L004_R2_001.fastq.gz

MM-265-DNA-T-01-01_L003_R1_001.fastq.gz
MM-265-DNA-T-01-01_L003_R2_001.fastq.gz
MM-265-DNA-T-01-01_L004_R1_001.fastq.gz
MM-265-DNA-T-01-01_L004_R2_001.fastq.gz

MM-465-DNA-N-01-01_L003_R1_001.fastq.gz
MM-465-DNA-N-01-01_L003_R2_001.fastq.gz
MM-465-DNA-N-01-01_L004_R1_001.fastq.gz
MM-465-DNA-N-01-01_L004_R2_001.fastq.gz

MM-465-DNA-T-01-01_L003_R1_001.fastq.gz
MM-465-DNA-T-01-01_L003_R2_001.fastq.gz
MM-465-DNA-T-01-01_L004_R1_001.fastq.gz
MM-465-DNA-T-01-01_L004_R2_001.fastq.gz

```

The final output needs to be as:

```
MM-465_normal_R1.fastq MM-465_normal_R2.fastq
MM-465_tumor_R1.fastq MM-465_tumor_R2.fastq
MM-265_normal_R1.fastq MM-265_normal_R2.fastq
MM-265_tumor_R1.fastq MM-265_tumor_R2.fastq

```

I've code below broken in three parts:

1. regular expression
2. process that iterates over files
3. workflow that calls the above two code blocks

regex:

```
def get_sample_id( file ) {
def matcher = (file =~ /(MM-\d+)-DNA-N-.+/)
    if (matcher) {
        return matcher[0][1]
    } else {
        return ''
    }
}

```

process:

```
process file_print {

debug true

input:

tuple val(sample_id), path(datasetFile)
val(sample_status)
val(read_value)

output:
publishDir '/arion/learn_new_concepts/', mode:'copy'

tuple val (sample_id), path("${output_file_name}")

script:

output_file_name = "${sample_id}_${sample_status}_${read_value}.fastq"

for (g in datasetFile){
println(g)

"zcat ${g} >> ${output_file_name} "

}

}

```

Workflow:

```
workflow {

normal_reads=Channel.fromFilePairs("/sc/arion/MM-*-N-*_L00*_R{1,2}_*.fastq.gz")
//MM-2658-DNA-N-01-01_L004_R2_001.fastq.gz  //

normal_forward_reads=normal_reads
        .map{ file -> return tuple( get_sample_id(file), file[1][0] )}.groupTuple(sort:true )

file_print(normal_forward_reads,"normal","R1").view() /// similar for reverse [1][1], R2 .. for tumor pass reads, "tumor", "R1"/R2"

}

```

I have not put code for reverse, but I get error as:

```
Caused by:   Process `file_print (1)` terminated with an error exit status (127) Command executed:   null

Command exit status:   127

```

I tried many things such as `${ output_file_name }`, `{ output_file_name }`, `$output_file_name` but in vain.

I'd like to improve few things in code:

a) The iteration/loop over file in `datasetFile` I think something like `each` could help, but I need order to be maintained

b) Order in which loop of `datasetFile` is created. For instance, I'd like to keep lanes L003, L004, L005 in this order. Or, L001, L002, L003 order, I hope this is clear.

I cannot move ahead or able to sort this by myself.

I do not quite understand how and when to use `.map{}`

Edited: after output added tuple for val and `path`. I still get same error.

# Answer
I think you are on the right track using `groupTuple(sort: true)`. Another way to get what you want would be to instead parse the `fromFilePairs` grouping key. Unless you're looking to do more than just concatenate (and decompress) FASTQ files, I think all the process really requires is an output prefix (in addition to the FASTQ files of course), for example:

```
params.reads = "/sc/arion/MM-*_R{1,2}_001.fastq.gz"

process concat_reads {

    tag { prefix }

    input:
    tuple val(prefix), path(reads, stageAs: 'reads/*')

    output:
    tuple val(prefix), path("${prefix}.fastq")

    """
    zcat -f ${reads} > "${prefix}.fastq"
    """
}

```

```
workflow {

    Channel.fromFilePairs( params.reads ) \
        | branch { key, fastqs ->

            def matcher = key =~ /^(MM-\w+)-DNA-(\w+)-.*_L(\d+).*$/

            matches: matcher
                def (_, sample_id, sample_type, lane) = matcher[0].findAll()

                def new_sample_id = "${sample_id}_${sample_type}"

                return tuple( new_sample_id, lane as int, *fastqs )

            other: true
                return tuple( key, *fastqs )
        } \
        | set { samples }

    samples.matches \
        | groupTuple( sort: true ) \
        | flatMap { sample_id, lanes, fq1s, fq2s ->
            [
                tuple( "${sample_id}_R1", fq1s ),
                tuple( "${sample_id}_R2", fq2s ),
            ]
        } \
        | concat_reads \
        | branch { prefix, fastq ->

            def idx = prefix.lastIndexOf('_')

            def sample_id = prefix.substring(0, idx)
            def suffix = prefix.substring(idx + 1)

            read1: suffix == 'R1'
                return tuple( sample_id, fastq )
            read2: suffix == 'R2'
                return tuple( sample_id, fastq )
        } \
        | set { result }

    result.read1 \
        | join( result.read2 ) \
        | map { sample_id, fastq1, fastq2 ->

            tuple( sample_id, [ fastq1, fastq2 ] )
        } \
        | view()
}

```

Results:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 23.04.1
Launching `main.nf` [jolly_noether] DSL2 - revision: 4a81f7a875
executor >  local (8)
[b8/a0b534] process > concat_reads (MM-265_T_R2) [100%] 8 of 8 ✔
[MM-465_T, [/path/to/work/9c/8b6fd328c848379423a0aba098cb74/MM-465_T_R1.fastq, /path/to/work/34/913780b1da902b5035f3bf56aac038/MM-465_T_R2.fastq]]
[MM-465_N, [/path/to/work/69/595120a1ee0a42217d06855ca9d1bc/MM-465_N_R1.fastq, /path/to/work/8f/31c406c5180073de577502fdb88fbf/MM-465_N_R2.fastq]]
[MM-265_N, [/path/to/work/1b/9c17cefb5aa51b3ee69832dc65d915/MM-265_N_R1.fastq, /path/to/work/d0/aac9f46e4a7f7165693b2537ad5630/MM-265_N_R2.fastq]]
[MM-265_T, [/path/to/work/d5/6295ff109df3ccb829a6f0a97b6e45/MM-265_T_R1.fastq, /path/to/work/b8/a0b53499ca13a523d811ecf772449b/MM-265_T_R2.fastq]]

```

Note that unless the FASTQ files are small, I would try to avoid decompressing the FASTQ files if possible. Most bioinformatics tools that deal with large FASTQ files will likely have gzip support, meaning you can replace the `zcat -f` with a `cat` in the example above (and remember to add the `.gz` extension to the output filename).
> 2 votes

---
Tags: nextflow
---


# Question
Title: GEO Methylation dataset "GSE73303" shows 0 features
I'm working for a project using GEO, in order to analyze some gene expression data. The SuperSeries I am referring to is GSE73304, but I found some problems working on the SubSeries GSE73303.

I started my analysis Microarray expression dataset and everything was fine. Now I have to investigate the other dataset, the one target as Bisulfite-Seq. But when I try to see inside it, I find there are correctly 20 patients, but 0 features, so seems it is empty. I attach R output:

```
library(GEOquery)
gsemet <- getGEO("GSE73303")

> gsemet
$GSE73303_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 0 features, 20 samples
element names: exprs
protocolData: none
phenoData
sampleNames: GSM1890515 GSM1890516 ... GSM1890534 (20 total)
varLabels: title geo_accession ... Sex:ch1 (44 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 27437034
28883895
Annotation: GPL11154

```

Could someone help me to find where I'm wrong? Or maybe, if the file is empty, pass me the right file?

Thank you very much in advance

# Answer
> 1 votes
It's NGS data, while GEOquery only returns data for arrays. You will either need to use the data they provude in the supplement of the GEO submission (if they do, it's not mandatory), or download and process fastq from scratch.

It's not gene expression by the way, it's methylation data, bisulfite sequencing.

---
Tags: r, bioconductor, genetics, geo
---


# Question
Title: Applying glmnet to identify predictors for subtypes
I am new to glmnet and other ML techniques, so apologise in advance if this sounds a trivial question. However, your guidance is very much appreciated.

So, I have nearly 700 genes with their expression profile across 4 conditions. I want to identify signature genes for each of the 4 conditions. I am using glmnet but I am not getting anywhere near the results.

so my data looks like this

```
           Subtype     HES4   FNDC10    PRKCZ PRKCZ.AS1
Sample1     NRT 5.452955 6.241387 8.546059  6.260245
Sample2     RT 5.125350 7.038719 8.900024  6.267325
Sample3     NRT 5.689560 5.903811 8.029434  4.213996
Sample4     NRP 5.834468 5.722589 8.061339  5.072274
Sample5     RP 5.401623 7.272111 9.108342  6.599629

```

purely in terms of numbers, the subtypes have following number of samples.

```
NRT 52
RT 42
NRP 30
RP 7

```

Q1. May be the groups are imbalanced and caution should be taken?

I have applied glmnet as following

```
modellinglasso=read.table("LRMatrix.txt", header=T, row.names = 1, sep="\t" )
modellinglasso$Subtype=as.factor(modellinglasso$Subtype)
modellinglasso$Subtype
#build a single pass (single fold) lasso-penalised model
lassoModel <- glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  family="multinomial")
plot(lassoModel, xvar="lambda")
Warning message:
In lognet(xd, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :
  one multinomial or binomial class has fewer than 8  observations; dangerous ground

```

```
#perform 10-fold cross validation
cv.lassoModel <- cv.glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  nfolds=10,
  family="multinomial",
  parallel=TRUE)
# plot variable deviances vs. shrinkage parameter, λ (lambda)
plot(cv.lassoModel)

```

```
idealLambda <- cv.lassoModel$lambda.min
idealLambda1se <- cv.lassoModel$lambda.1se
print(idealLambda); print(idealLambda1se)

[1] 0.1186125
[1] 0.1186125

# derive coefficients for each gene
co <- coef(cv.lassoModel, s=idealLambda, exact=TRUE)
head(co)

co.se <- coef(cv.lassoModel, s=idealLambda1se, exact=TRUE)
head(co.se)

# identify predictors for each sub-type
rownames(co.se$PR)[which(co.se$NRT != 0)]

[1] "(Intercept)"

```

output for

```
co.se$RT!=0
6 x 1 sparse Matrix of class "lgCMatrix"
            1
(Intercept) |
HES4        .
FNDC10      .
PRKCZ       .
PRKCZ.AS1   .
MMEL1       .

```

It doesn't matter which subtype I put, I get the same output "(Intercept)" instead of gene names.

Any help to get the gene signatures for the subtypes is sincerely appreciated.

Thank you

# Answer
Final answer. Keep in mind I work in Python so I'm trying to translate here.

1. The output of `coef(cv.lassoModel)` ... those are your genes of interest thats your **feature selection**. `head` I assume gives you the top values ... I think thats right (I need to check its not `coef(lassoModel)`) ... Basically thats where your results are.

This might be over simplistic but `co.se` without the `head` is your answer. Not every gene will carry a weight BTW. If the entirety of that is (700) `.` ... try point 2 below and then point 4. I don't think that will happen (see **Whats happened**)

2. I would redo this with `alpha=0` (ridge regression) and `alpha=0.5` elasticnet, but you don't need to
3. Generate an accuracy score

```
    y_predicted <- predict(best_model, s = best_lambda, newx = x)
    sst <- sum((y - mean(y))^2)
    sse <- sum((y_predicted - y)^2)
    rsq <- 1 - sse/sst
    rsq

```

Hopefully this is large e.g. around 0.8 (in Python this automatic BTW). You could ignore that, but its what I'd do

4. The intercept is simply the weight when your gene expression score is 0.
5. Lasso regression handles sparse data - so no worries about `"."`. I don't whether that Python's version of "nan" (not a number), or if its zero. ... @haci knows what this means and his answer is authoritative for the second part of your question.
6. In Python transformations are common, but that doesn't seem an option here, so find providing your accuracy score is reasonable, personally I'm happy.

In Python what you want is feature selection `print(pd.DataFrame({'Features': X_train.columns, 'Coefficients': lasso.coef_}))`

**Whats happened?** The coefficients for `head(co.se)` have been shrunk to zero, thats only the top 10? results, you've 700 weights. The full list *should* (famous last words) contain a minority of really big coefficients ... those are really important and those are your answer. So by printing only the top few results ... not the whole lot the results are still there you've just not seen them. Thats it, basically you were almost there and you could happily settle on that as the final answer. Basically, alot of genes ain't doing much a few genes are doing all the work - thats the answer. Below I stress repeatedly I don't think anything is wrong here and I personally like the result.

---

Background

**Lasso regression**, it's a cousin of ridge regression and key parameter is lambda (who'd guessed). What your plots are doing is tuning the parameter and finding the associated weight individually or collectively. Fairly straight forward. So like whats the value of lambda for your data set ('answered' below)? It's worth doing some reading on how lambda is calculated to understand the plot, if you're not sure how it works.

**Imbalance**, well to be honest I'm not convinced there's anything wrong with your analysis. Okay there's a lack of stationarity for some of your weights ... but some are fine and gist is pretty clear. Yeah okay there's a problem with stabilisation of lambda - thats not necessarily imbalance.

**Imbalance** is when it's 1 in 10 and over. You've got 1 in 6, 1 in 7, not huge. However, - *you might be right* \- and if you're worried simply repeat the analysis excluding `RP` (value 7). Then check whether the parameterisation plot (lambda against weight) is more stable or more "homogeneous" (multinomial plot), i.e. the genes behaviour in a more uniform way. If that is does work - it gives you nicer results AND `RP` is essential to the analysis thats a separate question about augmentation. That is complicated. For example, the multinomial plot with a 3 target categorical model might give a clear value for lambda.

**Transformation** Here, I would suggest if you want to stabilise the multinomial plot (if it's possible) I would look at a transformation. The data you presented doesn't look like it needs a transformation, but we can't see the whole data. In Python there's a *max-min* transformation - thats where I might start.

**Overall good results** Personally, I'd say these were good results, except the weights ain't homogeneous for some genes (most are fine BTW), but the patterns clear. The *red* (top), *blue* (bottom) and cyan (negative y-axis) are real results showing that within the 700 genes **3 genes** are having a massive impact. Sure you have positive and negative and thats how it works one gene is positively contributing to the classification, others are negatively contributing to it. Then you've the black and 'deep pink' lines which have stabilised around 5E4. **From 700 genes down to 5 genes**. Whats wrong with that?

**Whats the answer to the data set?**

Its ...

> NRT RT NRP RP

... is due to 'red', 'blue', 'cyan', 'black' and 'deep pink' genes. Thats it.

**Frank assessment** To be perfectly honest, looks good. No-one is going to say "Wooaaa lambda ain't stable, destroy! ;-)", they'll say "Wow thats cool and clear".

---

**How I'd do it..** I would look at basic ML robustness statistics and use this to assess the model rather than the parameterisation of lambda. Its horses for courses as they say. Upvotes.

**Note** I use Python not R for ML, both work well however. `glmnet` within ML is new to me BTW.
> 2 votes

---
Tags: r, rna-seq, machine-learning
---


# Question
Title: Read methylation data from tar files on GEO (GSE73303)
I'm trying to analyise methylation data using GEO, in order to conduct a multi-omics analysis. The SubSeries I am referring to is GSE73303. Unfortunately, the code used to inspect it produces empty file, so I downloaded tar files, as I found on bioconductor forum (ref. this).

I was able to download the tar files and extract thenm but afterwards I got a list of .bigwig files. I tried to read them on R, but actually I feel very lost and I can't understand how to trace back methylation data from these files.

If someone can guide me into this process, it would be really appreciated

Thanks in advance

# Answer
> 1 votes
As I've already told you here GEO Methylation dataset "GSE73303" shows 0 features GEO/GEOquery does not provide ready-to-use processed NGS data. It can be that the authors provide processed data, but if not (apparently not the case in this submission) then you have to start from the fastq files. Please see the many tutorials online on downloading fastq and processing bisulfite data, or try to use existing pipelines (nf-core methylseq might be an option) or seek a collaboration.

---
Tags: r, bioconductor, gene-expression, methylation, geo
---