text:
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?

Please format your error message/console output properly.

@RamRS I added full log file. Please If you want to help me with formatting my question, be more direct.

I was as direct as can be. What do you think I should do differently?

@RamRS I apologise for my comment. I did not mean to offend you. I was only seriously frustrated and I felt hopeless. This problem took me hours. I do not know, where I made mistake.

I'm not offended, Lukas. I was just wondering what I could have done differently.

Also, don't open files with notepad or Excel, try unix utilities or if you're on Windows, try notepad++

@RamRS Sorry, but I am new contributor. You can tell mi dircetly that I could my log message post fully and edit it little more. But still this is just pointless semantics.

@RamRS Thank you.

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

Thank you. Now I am really embaresed, that I opened those files in that way.


text:
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.

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.

source: Stack Exchange

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

sort: votes


text:
STACKS ref_map.pl no BAM records
I have 133 aligned and sorted BAM files that I am trying to run through the ref_map.pl script from the STACKS pipeline. My script looks as follows:
module load stacks

ref_map.pl -T 15 --popmap B2_filenames.txt --samples /home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/ --out-path /home/heather/scratch/ref_map_results/full_data

echo "ref_map.pl for B2 dataset completed"

I keep running into the issue where gstacks is unable to read the BAM files. This is what the log file looks like:
  Configuration for this run:
  Input mode: reference-based
  Population map: 'B2_filenames.txt'
  Input files: 133, e.g. '/home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/1401-47750.sort_filt.bam'
  Output to: '/home/heather/scratch/ref_map_results/full_data/'
  Model: marukilow (var_alpha: 0.01, gt_alpha: 0.05)

Reading BAM headers...
Processing all loci...
Error: No BAM records.
Error: (At the 0th record in file '/home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/1401-47754_B2.sort_filt.bam'.)
Aborted.

ref_map.pl: Aborted because the last command failed (1).
Last command executed by ref_map.pl was:
  /cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/stacks/2.64/bin/gstacks -I /home/heather/scratch/03_sort_and_map/filtered_and_sorted/B2/ -M B2_filenames.txt -O /home/heather/scratch/ref_map_results/full_data -t 15

I checked the BAM file with samtools - the headers look ok, and I don't get any output when I run samtools quickcheck, so I can't quite figure out what is causing the issue. I am following the same pipeline as I did with a previous batch of data and can't understand whats going wrong.

"the headers look ok, and I don't get any output when I run samtools quickcheck" - this would be true even if the files contained only headers, no? Do you see alignment info when you look beyond the headers?

I keep getting a similar error with the sorted and filtered bam file generated from samtools "Error: 'sample.bam': this is a SAM file (and BAM was specified).

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


text:
Protein-Protein Simulation using Gromacs
Can we do protein-protein simulation through Gromacs? if yes, can anyone teach me the correct procedure to get this done.
Thank you,

This is a very broad question. Could you kindly limit it to a specific technical issue? Secondly a reproducible example is always helpful. FYI the question has received a closed request.

I'd suggest looking at one of the many tutorials on the subject


text:
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

        """
        }

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.

Thank you. I removed def and it worked. My apologies I didn't quite understand: 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 What do I do? the example you shared is too difficult.

I have seen your reply few times on my posts as well on others, I don't quite follow queue channels, values channels link. What do I use when and for what? Sorry this is complicated then I can understand or/and follow.

@DeathMetal A better way, I think, would be to pass your sample_name variable (or  pair_id or whatever) down to your mutect2 process. It's likely you'll need this variable later in your workflow, so defining it in your workflow block likely makes more sense.

@DeathMetal I might be starting to sound like a broken record, but understanding how queue channels and value channels work should make writing workflows much easier. Queue channels can be consumed, but value channels can be read an unlimited number of times. In your example above, the intervals (targets) file and reference file should be passed in using value channels like in my previous answer. Note that value channels are implicitly created when a process is invoked with a simple value, like a file object.

Hi @Steve - can you please assist me on post https://bioinformatics.stackexchange.com/questions/21338/nextflow-merge-lane-files-from-samples-forward-and-reverse-read.


text:
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()

This might be better-suited to stackoverflow since it's primarily a selenium question rather than a bioinformatics question, but for what it's worth, I see the same behavior when I try your script.  (I've used selenium a bit before but never tried get for anything other than HTML.)  I don't suppose the e-utils support PDF paper download, do they?

To the best of my knowledge, E-utils primarily supports text, not PDF downloads. For some PMC-listed articles, the XML, mediline records lack full text, hence the need for browser automation. Though e-utils is my preferred choice, this situation demands an alternative. Nonetheless, the issue is now resolved

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)

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).

Thank you for your time! Unfortunately, it did not work, returning the error curl: (23) Failed writing received data to disk/application
Furthermore, I've noticed a potential issue. What if the article listed on PubMed is not available on PMC but on other platforms, such as the publisher's website? For example, PMID:37230009

You're welcome; just trying to help. That cURL error is very likely related to your curl install and file permissions (e.g. https://askubuntu.com/questions/1356327/cant-write-to-a-hidden-path-using-curl); you could try using wget instead, e.g. wget -O - https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=10026098 2> /dev/null | sed -n 's/.*\"\(ftp.*.pdf\)\".*/\1/p' should work. And as for articles on Pubmed that aren't on PMC, your question specifies 'articles on PMC', so that is what my answer addresses. If you want to download pdfs from any source, that is a totally different question.

Your solution is crashing because pdf="$(curl https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=$id_from_list | sed -n 's/."(ftp..pdf)".*/\1/p')" is empty for the pmcids that I used: 8968431, 8570779, 8290288, 8254498. But it worked for the pmcid that you provided. It seems pdfs corresponding to my pmcids are not present on the public-facing ftp server because when I tried your code with some of the pdfs that are listed there it worked. So, for my use case, browser automation is the only way so far

Ohh! Yep - I see the problem - only Open Access articles are accessible via the API (<OA><responseDate>2023-07-19 18:55:39</responseDate><request>https://www.ncbi.nlm.nih.gov/pmc/utils/oa/oa.fcgi?id=8968431</request><error code="idIsNotOpenAccess">identifier 'PMC8968431' is not Open Access</error></OA>). Thanks for taking the time to check!


text:
Bracken not working
I have some highly contaminated ancient DNA sequences. I have adapter removed and collapsed these and run them through Kraken2. The Kraken reports show multiple levels of taxa with good numbers of reads in each.
However, when I run Bracken - I keep getting an error saying it can't find any reads.
Error: no reads found. Please check your Kraken report

I have tried adjusting the taxonomic level (all the way up to Phylum) and get the same result. I have tried using different databases to run the Kraken analysis through and updating both packages.
Any ideas what could be going wrong?
Update:
Here's the head of one of the report files:

I am also using the latest version of Bracken downloaded via Bioconda. My command to run Bracken is via a bash script:
#SBATCH --partition=day
#SBATCH --output=slurm_bracken_job%J.out
#SBATCH --mem-per-cpu=10G
#SBATCH --cpus-per-task=32

# we load kraken2 and bracken into our environment:
ml kraken2
ml bracken

NAME=$1

srun bracken -d /workspaces/groups/database/nt-taxonomy-2021-02-04_braken/ -i ${NAME}.kreport -o ${NAME}.bracken -r 75 -l S

This is then run using:
for i in $(ls -1 ./*.report | sed 's/.report//'); do sbatch Test_Kraken.sh "$i"; done

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

It would be helpful if you could give details about your Bracken installation, the structure of your Kraken2 reports (what do the first few lines look like?), and the exact command you ran that threw the error

I've updated the main post - I have used a couple of different databases as was suggested in a github post that had the same problem but to no avail.

Thanks thats much better. My advice is one working first without the bash for loop. Bracken works, what I think is somethings not right in the run command.

I have just tried running a basic bracken -d /workspaces/groups/database/nt-taxonomy-2021-02-04_braken/ -i Sample_1.collapsed.report -o Sample_!.report.bracken -r 75 -l S and still the same error

I personally dunno. I do know is Bracken works and adjusts the frequency of the bacterial estimates from Kraken. I would personally try with their sample data which is here https://github.com/jenniferlu717/Bracken/commit/a6283e34c3925afe1e1bca3527e9fe0066fcae61#diff-223bc6965d2713e87fc5a41de1d0d94090ca3ea14ac92990c7191826b4e98df7

I had the same error. Here is how I fixed it. My problem was my kraken2 file. I ran kraken2 with the report zero counts option. When I deleted that, all was fine. I also used to run kraken2 with a confidence level of 0.6. I also removed that. Os you probably want to go back to your kraken2 report file and toy with how you generate this to make bracken work. Hope this helps.


text:
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?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

In your expectation, what is the DBA.2J strain tested against?

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.

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

I don't understand the difference between this and the vs offered by DESeq2 previously in your question. Given that the default test is Wald and betaPrior is TRUE by default for Wald tests, how is this result any different?

I'm not sure of the difference but this is the result I got when running the new code

Run it both ways and see if the number of genes or fold change value differs by much. I honestly don't think there should be any significant difference.


text:
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?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

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.

Thank you, the idea is using dask  distributed with multiple machines, then save parquet. But parquet files for a fastapi application is it good idea?

Please note, "accepts" (green tick-button) are much more welcome than thanks. Answered above


text:
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.

@M__ Thanks for the comment adjusted my question.

Added doing amplicon sequencing across different regions

Thanks thats slightly clearer. There is an issue here on the nomenclature. 15+K viral species is 6 times more viruses than designated by the ICTV (International Committee on the Taxonomy of Viruses), i.e. more viruses than the total number of species ever declared. Could you kindly define the term species? What you have stated is not taxonomically possible.

I see, using species not exactly clearly. They are all variants of the same virus. Will adjust

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.

Yes, they have controlled for some in vitro artefacts. I have plotted histograms of frequencies and each of the input variants. I got the expected picture you described, one variant taking a large portion of the reads a small number with a higher read count and a lot of low read count ones. I don't think I want to find a cutoff for the distribution, let's say I have one. I'm looking for a way to choose the number of input variants that maximizes observations over that threahold

Following up on your edit. For each input I have a control but I only have 4 inputs. It is too expensive/time consuming to sequence too many more. Do I just use the controls on the 4 inputs I have? Or is there some way I can simulate behavior for unknown inputs?

ok I think, I'm starting to pick up what you are saying. I have a few different options for what is considered a good observation. A hard threshold, item using some calculus (knee/inflection), as well as some methods based on spike ins (these are ideal probably).
But lets say for each input we can compute the fraction of observations above the threshold. Then fit a GLM to these 4 points? 4 seems a bit small for a GLM right?

Ok I think we are slightly talking about different things. I don't want to evaluate the controls or use them to select a threshold. I want to decide how many unique variants of my virus to sequence at once. Is that what your comment is describing and I'm just missing something?

Ok @TheNumber23 answered in comments above. Upvotes are encouraged.


text:
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?

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.

Thanks it worked.  I cannot understand why only one chromosome is being processed in lancet despite a channel being passed.


text:
Minfi and estimateCellCounts Error: in p[trainingProbes, ]: subscript out of bounds
I would like to plot using estimateCellCounts as I do not know my cell counts for the blood samples I am using however, it seem that I cannot get estimateCellCounts to work with EPIC arrays.
Detailing Attempts:
Attempt 1:

I first attempted it with EPIC arrays but it wouldn't take it.

estimateCellCounts(RGset,
                   compositeCellType = "Blood",
                   processMethod = "auto",
                   probeSelect = "auto",
                   cellTypes = c("Bcell",
                                 "CD4T",
                                 "CD8T",
                                 "Eos",
                                 "Gran",
                                 "Mono",
                                 "Neu",
                                 "NK"
                                ),
                   referencePlatform = "IlluminaHumanMethylationEPIC",
                   returnAll = TRUE,
                   meanPlot = TRUE,
                   verbose = TRUE)

And the associated error:
Loading required package: FlowSorted.Blood.EPIC
Loading required package: ExperimentHub
Loading required package: AnnotationHub
Loading required package: BiocFileCache
Loading required package: dbplyr

Attaching package: 'AnnotationHub'

The following object is masked from 'package:Biobase':
    cache

Warning message in data(list = referencePkg):
"data set 'FlowSorted.Blood.EPIC' not found"
Error in get(referencePkg): object 'FlowSorted.Blood.EPIC' not found
Traceback:
1. estimateCellCounts(RGset, compositeCellType = "Blood", processMethod = "auto",
 .     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
 .         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylationEPIC",
 .     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)
2. get(referencePkg)
3. get(referencePkg)

Adding the library("FlowSorted.Blood.EPIC") did not stop the error "data set 'FlowSorted.Blood.EPIC' not found".
Attempt 2:

Then, I looked up the error I got first on Github and found that converting to 450k was recommended. So I used the convertArray.
Then it threw up an error saying  "FlowSorted.Blood.450k" was not installed.

Attempt 3:
So I installed and loaded the library and re-ran to get the code below and the error below. However, this does not change the issue that it does not seem to run. The code I used and the error is below - any suggestions are appreciated:
library("FlowSorted.Blood.450k")

RGset_450k = convertArray(RGset,
                          outType = c("IlluminaHumanMethylation450k"),
                          verbose = TRUE
                         )

estimateCellCounts(RGset_450k,
                   compositeCellType = "Blood",
                   processMethod = "auto",
                   probeSelect = "auto",
                   cellTypes = c("Bcell",
                                 "CD4T",
                                 "CD8T",
                                 "Eos",
                                 "Gran",
                                 "Mono",
                                 "Neu",
                                 "NK"
                                ),
                   referencePlatform = "IlluminaHumanMethylation450k",
                   returnAll = TRUE,
                   meanPlot = TRUE,
                   verbose = TRUE)

Error:
Error in p[trainingProbes, ]: subscript out of bounds
Traceback:

1. estimateCellCounts(RGset_450k, compositeCellType = "Blood", processMethod = "auto",
 .     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
 .         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylation450k",
 .     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)

How would I make it get around this error as it does not output any results? Any suggestions/solutions are appreciated.

Seems that one person found that there was an issue with NAs ( https://support.bioconductor.org/p/9148370/ ) I am not sure if this is the case for me but adding this here as a comment just in case someone else has this issue with the NAs. The pull request seems to not have been approved yet.

I added source("/path/to/Forked_Minfi/R/estimateCellCounts.R") before the estimateCellCounts but I get the following error  after running the command ```Error in .isMatrixBackedOrStop(rgSet, "estimateCellCounts"): could not find function ".isMatrixBackedOrStop"
Traceback:

estimateCellCounts(RGset_450k, compositeCellType = "Blood", processMethod = "auto",
.     probeSelect = "auto", cellTypes = c("Bcell", "CD4T", "CD8T",
.         "Eos", "Gran", "Mono", "Neu", "NK"), referencePlatform = "IlluminaHumanMethylation450k",
.     returnAll = TRUE, meanPlot = TRUE, verbose = TRUE)```

I am the developer of FlowSorted.Blood.EPIC, I will try to provide you with some guidance.
First, the FlowSorted.Blood.EPIC is located on ExperimentHub, which means that you need to ask the program to download the database first before running estimateCellCounts (minfi) or estimateCellCounts2 (FlowSorted.Blood.EPIC).
For estimateCellCounts2, I explain this here if you are using a 450K array GitHub FlowSorted.Blood.EPIC. If you are using EPIC you can read the help file of estimateCellCounts2. Although we believe our method improves the cell estimates, you can load the library and use minfi's "auto" selection or the same option on estimateCellCounts2 (that is your decision as a user of the library).
library(FlowSorted.Blood.EPIC)
FlowSorted.Blood.EPIC <-
    libraryDataGet("FlowSorted.Blood.EPIC")
FlowSorted.Blood.EPIC

Then you can run estimateCellCounts according to minfi's help.
Good luck.

Welcome to the site @LucasSalas and thank you for your comprehensive response, much appreciated.

Thank you so much @Lucas Salas


text:
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

Thanks Sabrina, there is a recent Limma question which might help. Please could you supply the code you have tried and resulting error? I don't use Limma but it will almost certainly have a test data set. For example "16 Two-Color Case Studies" in the manual http://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf

Thank you so much for your help. I am pretty sure its the kind of data I am trying to input, I tried to use exprs() for my data and even that pulls up a data type error

Showing your code will help. There isn't many here who use Limma, but let's see. Also keep in mind not to confuse Limma with Limma-voom.

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.

Thank you very much for your help, I managed to find some data on GEO that I could practice on and its working with the code that I knew worked already since I had used it in a code along tutorial


text:
PLINK linear association testing command with an interaction and 2 covariates term
I am trying to write the right command for my association testing using PLINK 1.9 but I could not find a clear answer so far and not sure how to wrtie it correct.
I have genotype variants called and in PLINK bed formatted files under geno_data folder.
I have phenotype.txt file. First column as "Sample_ID", second column as "Treatment" (1 for Placebo and 2 for CompoundX) and third column has the Response variable as "ResX", which has integer number between 0 to 9.
I also have PCs file as cov.txt. The first and second columns are the same and includes Sample_IDs. The column names are "FID" and "IID". The third, fouth and fifth columns have PCs,with the column names "PC1", "PC2" and "PC3".
The Plink documentation did not help me as I am new in this. Therefore I need some examples command to clarify for various models as described below.
I need the correct PLINK (1.9) command (if not possible PLINK2 is also OK) for the following linear models:
Model 1:
ResX ~ Variant * Treatment + PC1 + PC2 + PC3

Second command I need is for this model:
Model 2:
ResX ~ Variant + Treatment + PC1 + PC2 + PC3

Model 3:
ResX ~ Variant * Treatment

Note I had asked this question on Biostar 11 days ago but did not receive any answer or comment yet.
Thank you


text:
Can files with different R1 and R2 lengths be trusted?
I received paired end amplicon sequence from a LAB with different lengths for R1 (320) and R2(280). Should I trust this lab to sequence other samples?
Also, I had to do a trimming over the R2 at 220(because Phred score was below 30) in the qiime2/DADA2 pipeline. Do you think trimming almost 60 bp from R2 won't result in false positives? I'm looking for precise results from my analysis. The lab claims this will give quality results as this is a specific experimental design of their lab. Should I trust the results I get from this data, also the lab? I have more samples to get sequenced.

Your title is not part of the question text; it's supposed to be a summary/primer. I'll fix it now but please do a better job next time.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

For paired-end Amplicon sequencing the lab has provided reads with different lengths, To use it for my downstream analysis I had to trim the R2 again, making it even shorter again. I have read that trimming too much of base pairs can result in false positives, during mapping of reads to 16s Database.  So I have 2 questions, first, the above statement of false positives is right or not, Second Do you think I can use this lab for further samples as well?

There's no reason why R1 and R2 must be the same length, what did you request from the sequencing lab? Also, plenty of labs use shorter sequences for their analyses. Without knowing more about what you are trying to use these sequences for, it is impossible to answer your question.

R1 and R2 can be different for different reasons. In your case looks like they did use different cycle of sequencing for R1 and R2. Har dot answer to your other question without knowing what the insert is or what the insert size is.

@Scot, I'm analyzing a person's gut microbiome. Read1 and Read2 are from the same paired-end amplicon sequencing run and originate from the same sample. Overlapping regions between Read1 and Read2 enable precise alignment and assembly, aiding identification of rare microbial species, detecting chimeric sequences, and improving alpha and beta diversity matrices. However, since R2 is shorter, I'm concerned about potential false positives in the final results due to limited overlap. Can you please advise? I'm new to Bioinformatics.

@Supertech Yes you are right, the Ct value is low for R1 and high for R2. So I guess concentration of R1 was higher than R2. Does the Ct value should remain consistent if  reads are from the same samples? Or Is it because of faults in experimental designs or condamination?

Hi @abi - I don't think anyone can tell you whether to trust your data or not - the key will be in how your analysis looks. But I can say that the trimming and the different lengths are not red flags to me. Much of the metagenomics work I've done used single end Illumina reads of 140nt, for reference.

Thanks, @Scot That's was helpful


text:
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:

regular expression
process that iterates over files
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.

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).

Thank you for your reply and code. I'm sorry, I cannot understand what you did in the code as it is scary and complicated for me. May I know what is wrong in my code? I've .fastq extension as a suffix. May be a little explanation of the steps you did would be highly helpful to me. For example, when to branch branch.  Or, how do I put normal in the concatenated files?

what did you do after concat_read  \ | branch what is this doing? Why do we need this?

Hi I've edited my code for output by providing tuple val and path in parenthesis. However I get same error. The error messages don't convey any meaning, they are incredibly cryptic. :(

@DeathMetal I suspect your error message is being truncated because it tells us exactly why it failed: .command.sh: line 2: null: command not found. Try running your workflow instead using nextflow run main.nf -ansi-log false. Since your script block returns null, this becomes the command Nextflow tries to run for better or worse. I think you could replace this with a script block containing: zcat -f ${datasetFile} > "${output_file_name}" to get it running. Also, use normal_reads.map{ key, files -> to unpack the fromFilePairs output.

@DeathMetal Note that processes do not necessarily process items in the order that they are received. For example, concat_reads can return results in any order. So, if we want to get back the reads where the R1 and R2 files are the first and second in the pair, respectively, then we need some way to separate the outputs so they can be joined correctly. In the above, I just used branch to create two channels: read1 and read2 and then a join to put them back together.

@DeathMetal Also, if you wanted to replace N with normal and T with tumor, I think the place to do this would when defining the new_sample_id in the above. You could use a Map for example.

Let us continue this discussion in chat.

I am unable to understand def (_, sample_id, sample_type, lane) = matcher[0].findAll() why there is a _ Is there a way to print/echo these outputs?

I want output concatenated from  concat_reads process only. They are generated in the folder as required. Is it necessary to write/do code below it where you join?

@DeathMetal If you just need the output of concat_reads, the code after and including the last branch statement is not necessary. But usually you would want to do something with these new pairs of files, which is why I included some code to show how this could be done. In reality, there are lots of ways to do this though. The underscore is just a placeholder for the full match. In Python and Groovy multi-assignment, at least, it is convention is to use an underscore for un-used variables. Use samples.matches.view() or insert a println(your_variable) into the branch statement to debug


text:
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

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.


text:
David tool online not converting most genes
I'm using DAVID tool online, in order to perform functional annotation chart. I extracted the most relevant genes from a dataset (GSE73304), but have the problem that DAVID is not able to understand the most part.
I tried to put my 250 most relevant genes, such as
A_33_P3230339
A_21_P0014879
A_23_P373568
A_23_P206371

into David, using as identifier "AGILENT_ID", as I found in the website. However, David is able to recognize just 89 of them (which are classified correctly as "homo sapiens"). I have the same problem if I try to do the gene_id conversion.
This might be not a David's problem, since I tried also on R withoud good results.
If someone knows the reason of this issue and how to solve it, it would be really appreciated
Thanks in advance

You may find annotation here : https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL17077

First try converting your Agilent probe IDs to Entrez ID or Gene Symbols and then run functional annotation analysis. You also need to understand that Gene Symbols or Entrez IDs are not available for all probe IDs.


text:
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

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

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)

I would redo this with alpha=0 (ridge regression) and alpha=0.5 elasticnet, but you don't need to
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

The intercept is simply the weight when your gene expression score is 0.
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.
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.

Thank you for your reply and detailed explanation. What I do not understand is why I am getting output as Intercept instead of gene names.

Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :(

Technically you are getting a sparse Matrix with co.se$RT!=0 and the "dots" are representing zeros or in this specific case FALSE values. You can see yourself by as.matrix(co.se$RT!=0).

All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples.

@Angelo I've addressed the comments in the first part of the answer.


text:
Filtering cohorts in PySpark returns Py4JJavaError
I'm working with UK Biobank data. I made a cohort with the Cohort Browser and tried to extract the data with this code  df = participant.retrieve_fields(names=field_names, filter_sql = cohort.sql, coding_values = "replace", engine=dxdata.connect()).
This is part of the text of the error: Py4JJavaError: An error occurred while calling 04100.showString. org.apache.spark. SparkException: Exception thrown in awaitResult:. The full error message is cut off.
However, if I remove the filter - df = participant.retrieve_fields(names=field_names, coding_values = "replace", engine=dxdata.connect()) – there is no issue.
What kind of error is this and why would filtering for a cohort help?

Thanks for the question and the edit. This is Python API for a Biobank - looks very interesting. Just to clarify the implication is you can connect to Biobank ok, but when you add a filter this stops the operation. However the code that connects okay has not been shown. Is that all correct?

Yes, that's correct

Okay, it's a bit hard without knowing the underlying code.

It looks like the issue comes from one of the fields I was filtering for. When I remove that filter the code works again. Still need to investigate why that happened.

Could you kindly expand on the answer please, specifically which field?

It really isn't clear. Possibly something in self-reported data, field 20002. I ended up using a workaround by filtering for fields separately and joining them with df.union(df2)

Generally remote biological servers, for example UCSC genome server, will not perform SQL queries for a free account. They'll allow downloads, but not queries. That is probably the same reason here.
The reason is simply these are often huge databases and queries can take up a lot of their computational resources, particularly if everyone is doing them. I don't know BioBank UK, but suggest its common practice.
If thats true the only difference with Biobank and UCSC is the UCSC server refusals are stated very clearly, "this operation ain't permitted for free accounts" (paraphrasing). What might be happening in an API error is not being picked up and relayed to you. Depends how the API is written so it's just reporting a timeout.

The account is paid, not free

@eclarette If you login directly you will receive a more informative message than going through an API, i.e. perform the SQL query there. The API isn't reporting the error message.


text:
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

If I've not already said this welcome to the site. The answer below is reasonable because generally we rarely provide overviews or manuals in the answer. What you need to do is break the question down to individual parts that don't work - divide and conquer - and provide the code for each stage. At present this isn't a cool question for this site because it's too broad. I hope that helps. In other words submit multiple focused questions rather than a very general question requesting an unspecified amount of help.

Ok I modify it based on your suggestion, thanks

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.

Yes I saw your reply thanks. I was just wondering, since I'm not an expert, if it was possilbe to follow the approach suggested in the reference that I linked.

@Riesz98 it would be etiquette to upvote.


text:
Error while opening xmfa generated by progressiveMauve
I used commandline progressiveMauve to align two genome .gbk files. The alignment was successful. But when I try to visualize the final .xmfa file with Mauve, I get the following error.
Exception in thread "Thread-1" java.lang.RuntimeException: Similarity index is too large.
at org.gel.mauve.histogram.ZoomHistogram.allocateIndex(Unknown Source)
at org.gel.mauve.histogram.ZoomHistogram.<init>(Unknown Source)
at org.gel.mauve.SimilarityIndex.<init>(Unknown Source)
at org.gel.mauve.XmfaViewerModel.init(Unknown Source)
at org.gel.mauve.XmfaViewerModel.<init>(Unknown Source)
at org.gel.mauve.ModelBuilder.buildModel(Unknown Source)
at org.gel.mauve.gui.FrameLoader.loadFile(Unknown Source)
at org.gel.mauve.gui.FrameLoader.run(Unknown Source)
at java.lang.Thread.run(Thread.java:745)

How can I solve this?
Is there another way to visualize the xmfa file.
My aim is to visualize an inversion in one of the genomes.


text:
How to subsample BAM based on read line number?
Consider a sorted, indexed, only mapped reads containing BAM file.
Is there a way to get a sub BAM based on read line numbers?
For example get a sub BAM file which contains read line numbers [1, 3, 8, 70] of the original BAM.
This can be done iteratively using a counter but its too slow. Is there a way to do this using samtools view or some tool like it.
Edit:

Size of the BAM file is ~7GB. Its human exome.

I'm using Julia (https://github.com/BioJulia/XAM.jl/) to iterate over the BAM and it takes me ~300 seconds to iterate over all the records. I feel this is too slow.

I'm running a function on each read of the BAM file and it decides which bucket the read has to go into. So this is specific subsampling.

I think samtools view - and awk route seems viable.

Welcome to the site . It would be helpful to state the size of the bam file, and particularly how you did this "reiteratively" (code if possible) - is essential because samtools is not the only solution. Like did you use Python and some sort of re, which would explain it? Also just run time, what is slow? Finally any additional information, e.g. it's human always helps, e.g. if it turns out to be bacteria ... thats different.

Please elaborate the purpose. How would you know that something in line x,y,z, is interesting so you want to subset for it? What is the context? Is it random subsampling or specific subsampling? In the simplest case you print the bam as sam via samtools view -h, skip header lines, and then make sed or awk print only lines overlapping an array of numbers you define upfront. That's almost no bottleneck, it just has to go through the entire file once which is not slow.

I've added more information in a edit. Can you please give me a sample code on how I can use AWK to print only lines overlapping an array of numbers.  I know to print Nth line you could do awk 'NR==N'

Julia is supposed to be fast (never used it). samtools will be faster.

Samtools+awk:
samtools view -h aln.bam \
  | awk 'BEGIN{while((getline<"lineno.txt")>0)a[$1]=1}/^@/{++n;print}!/^@/&&a[NR-n]' \
  | samtools view -b --threads=4 > sub.bam

I would recommend mawk which is much faster than other awk.
I believe using Julia should be faster as it would not need to convert BAM to SAM and then convert back. Five min to read through a 7GB BAM file sounds about right. That is not slow.


text:
bcftool merge of 5K samples makes ID a large string
I am using bcftools to merge 5000+ vcf files  using this command
bcftools merge -l list_of_vcfs.txt

I am getting the following error
 Warning: The row size is too big for FORMAT/PL at chr1:23899810, requires 12204615120 bytes, skipping.
The REF prefixes differ: C vs N (1,1)
Failed to merge alleles at chr1:53128382 in <sample>.vcf.gz

Any help on how to fix this issue would be appreciated.

Thanks for the edit @bioinfo_comp and welcome to the site. It's a good sign when the bot message is acted upon, particularly as quickly as you have done. Hope you stick around.


text:
Correctly coding indels for ped/map files in plink and VCF
I have a challenge with how to correctly code indels in plink (I plan to convert my binary file set to VCF and use it for lift over in downstream analysis).
For instance, my dataset has an indel with rsID rs72552763 and, alleles ATGAT and AT. Its position is chr6: 160139849 - 160139853 (GRCh 38.p14).
In the .map file, which bp position should I indicate?
With respect to the .ped file, how should I code the alleles/genotypes?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Welcome to the site. I would simply recommend referring to the manual https://www.cog-genomics.org/plink/ and https://zzz.bwh.harvard.edu/plink/dist/plink-doc-1.07.pdf . There is plenty of information on both .ped and .map therein to address format concerns, in addition there will be reproducible examples. It would be useful to specify the sort of command you wish to run e.g. plink --dosage myfile.dat --fam mydata.fam --map mymap.map.