text:
Pathway Enrichment Analysis for one specific pathway?
I have a gene set containing a set of 1500 genes in which I did differential gene expression using limma voom pipeline. After this step I wanted to analyze a specific set of pathways the metabolism set from KEGG database.
hs_kegg_df <- msigdbr(species = "Homo sapiens") %>%
 dplyr::filter(
gs_cat == "C2", # This is to filter only to the C2 curated gene sets
gs_subcat == "CP:KEGG" # This is because we only want KEGG pathways
)

Then since I only want to analyze metabolism I did this (keep only the pathways envolved in metabolism):
hs_kegg_df_with_gs_source <- subset(hs_kegg_df, grepl("^hsa00|^hsa01|^hsa02", gs_exact_source))

My problem is that, for example, for the pathway  KEGG arginine and proline metabolism the set contains 54 genes in which my kit only covers 4, and these for 4 genes are Differentially expressed yet when doing: ( I decided to do ORA because with 1500 genes GSEA is not recommended).
em <- enricher(genes, minGSSize=1, maxGSSize = 20,  universe = names(gene_list), TERM2GENE=hs_kegg_df,pAdjustMethod = "fdr")

This set does not appear with significant p value. I also tried passing to the enricher just that pathway in particular, the p value is 1. Because in hypergeometric test, the background genes are the same as the as DGE genes in all pathways giving this way 1-0=1.
How can I analyze just one pathway in particular given this specific conditions?
These are the 4 genes DGE and included in my kit:
  gs_name                              gene_symbol
  <chr>                                <chr>
1 KEGG_ARGININE_AND_PROLINE_METABOLISM ALDH2
2 KEGG_ARGININE_AND_PROLINE_METABOLISM ARG2
3 KEGG_ARGININE_AND_PROLINE_METABOLISM NOS1
4 KEGG_ARGININE_AND_PROLINE_METABOLISM ODC1

And this are all the 59.
here all the genes dge and all the genes in pathway:
 genes_dge <- c("ALDH2", "ARG2", "NOS1", "ODC1")
 all_genes <- c("ACY1", "AGMAT", "ALDH18A1", "ALDH1B1", "ALDH2", "ALDH3A2", "ALDH4A1", "ALDH7A1", "ALDH9A1", "AMD1", "AOC1", "ARG1", "ARG2", "ASL", "ASS1", "AZIN2", "CKB", "CKM", "CKMT1A", "CKMT1B",s"CKMT2", "CPS1", "DAO", "GAMT", "GATM", "GLS", "GLS2", "GLUD1", "GLUD2", "GLUD2", "GLUL", "GOT1", "GOT2", "LAP3", "MAOA", "MAOB", "NAGS", "NOS1", "NOS2", "NOS3", "OAT", "ODC1", "OTC", "P4HA1", "P4HA2", "P4HA3", "PRODH", "PRODH2", "PYCR1", "PYCR2", "PYCR3", "PYCR3", "SAT1","SAT2", "SMS", "SRM")

Thank you for the answer @M__  its samples from two different subtypes of human cancer. And its "C2".

Two questions

Significance
Lack of consistency between enzymes (log2) within the arginine and proline metabolism pathway.

.  1. We've dealt with significance in a past post on Limma (it was not voom though). I'll try and find the post and link it here.

The answer is straight forward: whilst all four enzymes include Arginine and proline metabolism pathway

only ARG and ODC1 share exactly the same pathways. Thus these should give very similar results.

NOS1 is a key enzyme in 6 different pathways.

ALDH2 is involved in 13 different pathways.

Therefore NOS1 and ALDH2 would have "confounders" (from an experimental perspective, not host metabolism) because they are more complex than singly arginine and proline metabolism pathways.
The other factor is in cancer patients the transcriptomics will get messed up, ALDH2 could go all over the place.
Prediction  ARG = ORD < NOS1 < ALDH2
**ALDH2** aldehyde dehydrogenase 2
Pathways
Glycolysis / Gluconeogenesis
Ascorbate and aldarate metabolism
Fatty acid degradation
Valine, leucine and isoleucine degradation
Lysine degradation
**Arginine and proline metabolism**
Histidine metabolism
Tryptophan metabolism
beta-Alanine metabolism
Glycerolipid metabolism
Pyruvate metabolism
Pantothenate and CoA biosynthesis

**ARG** arginase
Pathways
**Arginine and proline metabolism**
Glutathione metabolism

**NOS1** nitric oxide synthase
Pathways
Arginine biosynthesis
**Arginine and proline metabolism**
Calcium signaling pathway
Phagosome
Apelin signaling pathway

**ODC1** ornithine decarboxylase 1
Pathways
**Arginine and proline metabolism**
Glutathione metabolism


text:
Plot error and blank tab of plot
I am trying to plot by using these commands , but this is the error I am getting, and I tried multiple times, even by making the window size larger, the [lot tab is empty/blank too.
plot1 <- VariableFeaturePlot(nsclc.seurat.obj)
LabelPoints(plot = plot1, points = top10, repel = TRUE)***
*When using repel, set xnudge and ynudge to 0 for optimal results
Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  :
  Viewport has zero dimension(s)
In addition: Warning message:
Transformation introduced infinite values in continuous x-axis*

Welcome @Samia. Could you kindly clarify your study? I get the impression this is a Seurat analysis and I assume you are plotting the PCA. Any further information on the study would help understand the error output. Basically that the PCA is empty "zero dimension(s)" - not good. The error is likely further upstream.

source: Stack Exchange

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

sort: votes


text:
Paired end reads with different read lengths
I downloaded some raw RNA-seq reads from NCBI using the sra-toolkit (Accession numbers: SRR2053159-64). The reads were generated using ABI SOLiD, and in colour space. I downloaded them in base space in fastq format using the -B option in fastq-dump. But after that, when I am looking at the data after running FastQC, I find that the forward read length is 50 while the reverse reads length is 35. How is this possible? This is causing all the reverse reads to fail the filtering criteria set in Trimmomatic MINLEN:36. What is the correct way to work around this?
I have already set LEADING:20 TRAILING:20 AVGQUAL:25 SLIDINGWINDOW:10:30 MINLEN:36 Would it be a good idea to entirely remove the MINLEN:36, or would it be better if I make the threshold 25?
Actually I am trying to compare the transcriptomes of different bacteria under a certain stress using raw reads from NCBI. While for the rest of the set of bacteria the sequencing was done using either Illumina or Solexa, this one was done using ABI SOLiD. Is it a good idea to include this ABI SOLiD data along with the others for analysis?

I find that the forward read length is 50 while the reverse reads length is 35. How is this possible?

Forward cycle length was longer on SOLiD than reverse cycle length; this is common for SOLiD data. If you want to work with your data, you'll need to adjust your trimming settings appropriately to account for this (e.g. getting rid of the minimum length entirely, or setting it to a lower value). Trimmomatic does have other trimming options that are not length-based. Given the options that you've mentioned, it should be fine to just drop that minimum length parameter.

Is it a good idea to include this ABI SOLiD data along with the others for analysis?

No, especially not if you're converting to base space, then mapping. SOLiD data should be mapped in colorspace, so that small changes in base transitions don't have a huge effect on the mapped sequence. SOLiD data converted to base space sequence is not appropriate for mapping because it mixes different error modes (i.e. sequencing error, and base transition error). Bowtie (as in Bowtie v1) can map to a colorspace reference, and will correct errors based on that reference in the resultant SAM output. This means that if you trust your reference, you can use Bowtie as a data cleaner for colorspace data.
However, on a practical level, the output of ABI SOLiD machines for existing datasets was low enough that (coupled with the bioinformatics issues) it makes more financial sense to resequence samples than to try to navigate around the pitfalls of colorspace representation.
I have previously offered the following steps for dealing with colorspace data:

Transfer all the colorspace files onto an external hard disk
Delete all other copies of the colorspace files
Remove the hard drive from the computer
Use a sledgehammer or similar to squash the disk platters closer together
Withdraw \$500 from the bank
Place the \$500 on top of the hard drive
Return the hard drive (with the money) back to the client
Report to the client that there was insufficient data for a suitable analysis, and recommend that the experiment is repeated using a base-space sequencer

A More detailed Explanation
[from here]
In colorspace sequences that are converted to base space, any error in the sequence will cause all the following bases to be incorrect. Any colour-space to base-space conversion needs to take into account (and correct) errors so that the base-space sequences are correct. When there is a sequence difference, the conversion needs to make sure that only that position is changed in the base-space version.
Consider the following sequences that map to the same position:
reference: G101320112
sequence1: X101120112
sequence2: X101312011
sequence3: X201320312

[I chucked an INDEL in there to make things a bit harder]
A naive base-space conversion would convert these sequences as follows:
reference: GTTGCTTGTC
sequence1: GTTGTCCACT
sequence2: GTTGCAGGTG
sequence3: GAACGAATGA

[apologies if my conversion is incorrect. Fixes appreciated]
Very similar colour-space sequences, but very different base-space sequences.
A more correct conversion would notice where the errors were in the sequences relative to the index, and modify the next colour-space base as well to something that looks appropriate:
reference: G101320112
sequence1: X101100112
sequence2: X101313011
sequence3: X231320332

This would end up with these converted sequences:
reference: GTTGCTTGTC
sequence1: GTTGTTTGTC
sequence2: GTTGCATTGG
sequence3: GATGCTTATC

Which look considerably better.
I hate colorspace because the conversions are very unintuitive, and difficult to explain to other people. About the only nice thing is that reverse complement is just the reverse, but this also means that aligners and assemblers need to be modified to account for that when working in colorspace (or double-encoded colorspace), and you can get weird unexpected chimeras (e.g. poly-A tails and poly-T heads merging). You can save a lot of pain and confusion by sticking with a base-space sequencer.

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


text:
How to reverse complement a genbank file?
I am visualizing the alignment of a particular genomic region across different genomes with mauve. In certain genomes the region of interest is in the negative strand and in others the positive strand. What tool can I use to reverse complement my genbank files without losing the features and annotation?

I have not really used Mauve, so I don't know whether it might have an in-built setting for that.
SnapGene Viewer should be able to do that - it can open genbank files and export them, and View > Flip Sequence gets you the reverse complement without affecting the annotations. SnapGene can also do that; I only specified Viewer because it's free.

Thank you. It worked.

From the biology side of bioinformatics and the easier option, I would use Mauve.  If you want to make the project more specific and tailored to your needs, and you have python skills, I would write Python code to complete this task.  I would specifically use the bioinformatics Python library called Biopython: https://biopython.org
pip install biopython
The second option is more involved but it gives you more control on your project/research since you're the one creating it.  I personally would go with the 2nd option if you have coding skills.  If you don't have coding skills, stick to the 1st option.  Im sure there are also tools you can use that are found on Galaxy:
https://usegalaxy.org
I am not sure what tools there are that would be good for your task and needs, but from using Galaxy in the past, im sure you can find something! :)


text:
GTF To Use With Broad Institute Fasta In STAR
To make a STAR reference genome used for alignment you need a genome fasta (the resources fasta provided by the Broad Institute in my case) and a GTF fasta. I tried to use the Ensembl human genome fasta fro hg38 but that does not work. I have looked in the GATK resources package but there is no .gtf file there. I really don't know what to do in this case? Is there anyone who can provide some guidance as to what .gtf file I should use?

If you are looking for a genome and its corresponding annotations, you can use the following, it is from nf-core (link to the original doc):
#!/bin/bash

VERSION=108
wget -L ftp://ftp.ensembl.org/pub/release-$VERSION/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz
wget -L ftp://ftp.ensembl.org/pub/release-$VERSION/gtf/homo_sapiens/Homo_sapiens.GRCh38.$VERSION.gtf.gz

BTW the most up-to-date version as of today is 109.

When I use the Ensembl it doesn't work saying Fatal INPUT FILE error, no valid exon lines in the GTF file: and then the path to the Ensembl gtf.

That makes no sense, Ensembl's GTF is the most comprehensive GTF there is.

Thank you @gringer

I've had success in the past generating STAR indexes using the Gencode GTF files available from here:
https://www.gencodegenes.org/human/

Comprehensive gene annotation (ALL regions)
Genome sequence (ALL)

If you make sure that you're getting both the genome and the GTF files from the same source, it's more likely that they will work well together.
Note that the ALL annotation includes assembly patches and alternate haplotypes. The alternate haplotypes mean that population variation is captured, which usually increases the alignment rate of reads.
Unfortunately, most bioinformatics workflows don't properly handle alternate assemblies. As one example, reads aligning to those areas may end up classified as mapping to multiple genomic locations, even though they actually aren't. So it's a tradeoff of a more accurate alignment vs more hassle dealing with those alignments downstream.


text:
Analysis of RNA sequencing with zUMIs
I wrote the .yaml file which was according the instruction from zUMIs:
###########################################
#Welcome to zUMIs
#below, please fill the mandatory inputs
#We expect full paths for all files.
###########################################

#define a project name that will be used to name output files
project: SmartSeq3

#Sequencing File Inputs:
#For each input file, make one list object & define path and barcode ranges
#base definition vocabulary: BC(n) UMI(n) cDNA(n).
#Barcode range definition needs to account for all ranges.You can  give several comma-separated ranges for BC & UMI sequences, eg.BC(1-6,20-26)
#you can specify between 1 and 4 input files

sequence_files:
  file1:
    name: /home/isidora/IRCCS_Candiolo/Fastq/*_R1_001_val_1.fq ##/home/isidora/IRCCS_Candiolo/Fastq/zUMIsFiles1/ #path to first file
    base_definition:
      - cDNA(1-150)
      - BC(1-8)
      - UMI(12-19)
  file2:
    name: /home/isidora/IRCCS_Candiolo/Fastq/*_R2_001_val_2.fq ## /home/isidora/IRCCS_Candiolo/Fastq/zUMIsFiles2/ #path to second file
    base_definition:
      - cDNA(1-150)
      - BC(1-8)
      - UMI(12-19)

#reference genome setup
reference:
  STAR_index: /home/isidora/IRCCS_Candiolo/Fastq/ #path to STAR genome index
  GTF_file: /home/isidora/IRCCS_Candiolo/Fastq/hg38.refGene.gtf #path to gene annotation file in GTF format
  exon_extension: no #extend exons by a certain width?
  extension_length: 0 #number of bp to extend exons by
  scaffold_length_min: 0 #minimal scaffold/chromosome length to consider (0 = all)
  additional_files: /home/isidora/IRCCS_Candiolo/Fastq/hg38.fa #Optional parameter. It is possible to give additonal reference sequences here, eg ERCC.fa
  additional_STAR_params: pGe.sjdbOverhang > 0 #Optional parameter. you may add custom mapping parameters to STAR here
  #output directory
out_dir: /home/isidora/IRCCS_Candiolo/Fastq/zUMIs/ #specify the full path to the output directory

###########################################
#below, you may optionally change default parameters
###########################################

#number of processors to use
num_threads: 20
mem_limit: null #Memory limit in Gigabytes, null meaning unlimited RAM usage
#barcode & UMI filtering options
#number of bases under the base quality cutoff that should be filtered out.
#Phred score base-cutoff for quality control.
filter_cutoffs:
  BC_filter:
    num_bases: 3
    phred: 20
  UMI_filter:
    num_bases: 2
    phred: 20

#Options for Barcode handling
#You can give either number of top barcodes to use or give an annotation of cell barcodes
#If you leave both barcode_num and barcode_file empty, zUMIs will perform automatic cell barcode selection for you!
barcodes:
  barcode_num: null
  barcode_file: /home/isidora/IRCCS_Candiolo/Fastq/SampleSheet.txt
  barcode_sharing: null #Optional for combining several barcode sequences per cell (see github wiki)
  automatic: yes #Give yes/no to this option. If the cell barcodes should be detected automatically. If the barcode file is given in combination with automatic barcode detection, the list of giiven barcodes will be used as whitelist.
  BarcodeBinning: 1 #Hamming distance binning of close cell barcode sequences.
  nReadsperCell: 100 #Keep only the cell barcodes with atleast n number of reads.
  demultiplex: yes #produce per-cell demultiplexed bam files.

#Options related to counting of reads towards expression profiles
counting_opts:
  introns: yes #can be set to no for exon-only counting.
  intronProb: no #perform an estimation of how likely intronic reads are to be derived from mRNA by comparing to intergenic counts.
  downsampling: '0' #Number of reads to downsample to. This value can be a fixed number of reads (e.g. 10000) or a desired range (e.g  10000-20000) Barcodes with less than <d> will not be reported. 0 means adaptive downsampling. Default: 0
  strand: 0 #Is the library stranded? 0 = unstranded, 1 = positively stranded, 2 = negatively stranded
  Ham_Dist: 1 #Hamming distance collapsing of UMI sequences.
  velocyto: no #Would you like velocyto to do counting of intron-exon spanning reads
  primaryHit: yes #Do you want to count the primary Hits of multimapping reads towards gene expression levels?
  multi_overlap: no #Do you want to assign reads overlapping to multiple features?
  fraction_overlap: 0 #minimum required fraction of the read overlapping with the gene for read assignment to genes
  twoPass: yes #perform basic STAR twoPass mapping

#produce stats files and plots?
make_stats: yes

#Start zUMIs from stage. Possible TEXT(Filtering, Mapping, Counting,Summarising). Default: Filtering.
which_Stage: Filtering

#define dependencies program paths
samtools_exec: samtools #samtools executable
Rscript_exec: Rscript #Rscript executable
STAR_exec: STAR #STAR executable
pigz_exec: pigz #pigz executable

#below, fqfilter will add a read_layout flag defining SE or PE

When I start the program with this code zUMIs/zUMIs.sh -c -y ~/IRCCS_Candiolo/Fastq/zUMIs/zUMIs.yaml, I get errors:
Using miniconda environment for zUMIs!
 note: internal executables will be used instead of those specified in the YAML file!

 You provided these parameters:
 YAML file:     /home/isidora/IRCCS_Candiolo/Fastq/zUMIs/zUMIs.yaml
 zUMIs directory:               /home/isidora/IRCCS_Candiolo/Fastq/zUMIs
 STAR executable                STAR
 samtools executable            samtools
 pigz executable                pigz
 Rscript executable             Rscript
 RAM limit:   null
 zUMIs version 2.9.7e

Wed Aug 23 12:06:25 CEST 2023
WARNING: The STAR version used for mapping is 2.7.3a and the STAR index was created using the version 2.7.4a. This may lead to an error while mapping. If you encounter any errors at the mapping stage, please make sure to create the STAR index using STAR 2.7.3a.
Filtering...
sh: 1: Syntax error: Unterminated quoted string
sh: 1: Syntax error: Unterminated quoted string
Wed Aug 23 12:06:26 CEST 2023
Error in eval(bysub, parent.frame(), parent.frame()) :
  object 'XC' not found
Calls: cellBC -> [ -> [.data.table -> eval -> eval
In addition: Warning message:
In data.table::fread(bccount_file, header = FALSE, col.names = c("XC",  :
  File '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//SmartSeq3.BCstats.txt' has size 0. Returning a NULL data.table.
Execution halted
Mapping...
[1] "2023-08-23 12:06:28 CEST"
Warning message:
In data.table::fread(cmd = paste(samtools, "view", filtered_bams[1],  :
  File '/tmp/RtmpKSSqU2/file57f5393a8fde' has size 0. Returning a NULL data.table.

EXITING because of fatal PARAMETERS error: pGe.sjdbOverhang <=0 while junctions are inserted on the fly with --sjdbFileChrStartEnd or/and --sjdbGTFfile
SOLUTION: specify pGe.sjdbOverhang>0, ideally readmateLength-1
Aug 23 12:07:18 ...... FATAL ERROR, exiting
Wed Aug 23 12:07:18 CEST 2023
Counting...
[1] "2023-08-23 12:07:30 CEST"
Error in fread(paste0(opt$out_dir, "/zUMIs_output/", opt$project, "kept_barcodes_binned.txt")) :
  File '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//zUMIs_output/SmartSeq3kept_barcodes_binned.txt' does not exist or is non-readable. getwd()=='/home/isidora/IRCCS_Candiolo/Fastq/zUMIs'
Execution halted
Wed Aug 23 12:07:30 CEST 2023
Loading required package: yaml
Loading required package: Matrix
[1] "loomR found"
Error in gzfile(file, "rb") : cannot open the connection
Calls: rds_to_loom -> readRDS -> gzfile
In addition: Warning message:
In gzfile(file, "rb") :
  cannot open compressed file '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//zUMIs_output/expression/SmartSeq3.dgecounts.rds', probable reason 'No such file or directory'
Execution halted
Wed Aug 23 12:07:32 CEST 2023
Descriptive statistics...
[1] "I am loading useful packages for plotting..."
[1] "2023-08-23 12:07:33 CEST"
Error in data.table::fread(paste0(opt$out_dir, "/zUMIs_output/", opt$project,  :
  File '/home/isidora/IRCCS_Candiolo/Fastq/zUMIs//zUMIs_output/SmartSeq3kept_barcodes.txt' does not exist or is non-readable. getwd()=='/home/isidora/IRCCS_Candiolo/Fastq/zUMIs'
Execution halted

I am new at bioinformatic field and I do not understand the meaning of these errors and the things I should correct.
Can anyone help me please?

I don't see any obvious missing ' or " (that's what "Unterminated quoted string" means) in the zUmis.sh version on github, but maybe yours is different? Can you changing downsampling: '0' to downsampling: 0 in your yaml file and relaunching? I doubt this is the issue, but it's the only quote I see there and maybe the script isn't picking it up correctly (you normally don't need to quote integers, and this is the only quoted one).

Also, should STAR_index: be pointing to a file instead of a directory? And did you see the warning? WARNING: The STAR version used for mapping is 2.7.3a and the STAR index was created using the version 2.7.4a. This may lead to an error while mapping. If you encounter any errors at the mapping stage, please make sure to create the STAR index using STAR 2.7.3a. Have you tried that?

I'd advise you to keep away from ZUMIs unless you're doing the same workflow that it was designed for. I spent many days trying to sort out issues with that program, and eventually gave up; it wasn't created with the expectation of errors or unexpected input, and had too many programming minefields for the data I was working with.
ZUMIs should be stopping on the first error and explaining what's wrong. As demonstrated by your situation, it's not; it's trying to keep on processing, barreling through and dragging errors along for the ride. It's not exactly beginner-friendly code.
In the end, I changed to using Salmon/Alevin, which had a more fail-resistant interface, and (eventually) similar barcode flexibility:
https://combine-lab.github.io/alevin-fry-tutorials/2021/custom-geo/

0
Thank you for your quick answer. I did the SmartSeq 3 protocol and it is advised to do bioinformatic analysis with zUMIs pipeline because of the specificities of results. Do you know is the Salmon/Alevin suitable for this protocol?

I don't know, sorry. If you've been advised to use ZUMIs, then do that, and try to get more help from the person who gave you that advice.


text:
Model optimization causes overfitting
I'm working on an ML model to predict colorectal cancer using data from the human microbiome (collected from fecal samples).
My dataset is made up of metadata (containing phenotypic information on the N=156 population on which the study is being carried out,
 and an OTU table containing a taxonomy of bacteria found and their abundance for each sample from one person (each person has one and only one sample).
I filtered the OTU table by eliminating sparse data (I eliminated OTUs below a certain threshold) after which I calculated the relative abundance for each OTU.
I deleted samples with missing values.
I made the metadata join (Subject ID    Sample ID    Age (years)    Gender    BMI (kg/m²)    Country of Residence    Diagnosis    TNM Stage    AJCC Stage    Localization) and the OTU table.
In the end, I ended up with 114 rows × 221 columns .
I normalized the data and then built an SVM(SVC) model with the rbf kernel with C=0.05 and gamma=0.001 I got
Training Accuracy: 0.967032967032967
Test Accuracy: 0.782608695652174
SVC(AUC=0.64)

When I try to optimize the model using the code below, the values I get for gamma and C give me overfitting (trainset accuracy is 100%).
# import packages
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
# range of values to try for gamma and C
Cs = [0.00001,0.0001,0.001,0.1,1,2,5,10,100,1000]
gammas = [0.00001,0.0001,0.001,0.1,1,10,100,10000]
# number of segments for cross validation
nfolds = 5 # number of segments for class validation
param_grid = {'C':Cs ,'gamma': gammas}
# search for the best parameters using cross-validation
grid_search = GridSearchCV(SVC(kernel='rbf'),param_grid,cv=nfolds)
grid_search.fit(X_train,y_train)
grid_search.best_params_

It seems like an overfitting problem, a sample of the dataset would be appreciated, I think the model might be hanging onto a few features to make the prediction.

You are missing a validation split triggering data leakage. The error is here:
grid_search.fit(X_train,y_train)

You can't use X_train and y_train here. The correction is
grid_search.fit(X_val,y_val)

Generating X_val and y_val is described below.
The code used to generate

Training Accuracy: 0.967032967032967
Test Accuracy: 0.782608695652174

...  included one test_split_train and you need to perform a second test_split_train to obtain X_val, and y_val. The original code gave a pretty good as a result, the training score will always be higher.
The diagnosis of the code presented is there is 100% data leakage. What is missing is a second test_train_split. Again one test_train_split is present in the original code and there would have been a,
import pandas as pd
from sklearn.model_selection import train_test_split

X = # a dataframe containing the features remove the target y (below)
y = # a same dataframe but a single column which is the target
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=1)

The above will give a training proportion at 80% and test proportion at 20%. It's called an 80:20 split. That code will have been in the original code for the first calculation, I know that because the results are pretty standard for a routine ML test_train_split operation. Whats happening is the training set is being held absolutely separately to the test set. There is 100% compartmentalisation between the two data sets. Thats critical because if they got confused the calculation falls over, you will be testing the exact data thats just been trained. That confusion would be term "data leakage".
In the code you presented what has happened is the same data used for parameterisation is exactly the same data for training (X_train, y_train). Again it's called data leakage and parameter grids are very sensitive to this.
To correct this when a parameter grid is used it must have its own data set. That is instead of having a X_train,X_test,y_train,y_test the test_train_split is modified to generate X_train,X_test,y_train,y_test, X_val, y_val. The code to generate this is,
import pandas as pd
from sklearn.model_selection import train_test_split
 X_train, X_test, y_train, y_test
    = train_test_split(X, y, test_size=0.2, random_state=1)

 X_train, X_val, y_train, y_val
    = train_test_split(X_train, y_train, test_size=0.25, random_state=1)

This approach is singly leveraging test_train_split by splitting the training data set twice. This example creates a 60:20:20 split. Some say that 70:20:10 split is used for a parameter split (10 for the parameter search). Matter of opinion (you just adjust the proportions).
Thus

the grid search is strictly for the parameter split (X_val, y_val)
and those parameters at then using in training split (X_train, y_train)
the testing is then performed (X_test, y_test)

What this process does is ensure no part of X_val, y_val is ever seen in the training set. If that happens again - data leakage. Here only the optimised parameters are seen by the training set.

Summary Introducing validation parameter split into the calculation isn't trivial and there is nothing wrong with your original result. Thus you don't necessarily need it. Generally speaking 80% accuracy is fine, yours is 1% lower (can be overlooked).
Parameter tuning will increase the accuracy score but it will be 85% rather than 79% and if it's really gone well then it might be 90%. What I'm saying is it will improve things, but will not absolutely solve the accuracy issue, it's just a general and accepted issue within ML.

Just to finally add there are biological reasons in cancer prediction why accuracy is a key parameter, but thats a different question.

Comments: causes lower scores
What I'm saying the correct approach. I recognise that Geek for Geeks uses your approach - but they are wrong.
A separate validation split will certainly give result in lower scores than grid_search.fit(X_train,y_train) thats inevitable, but it is correct in stopping overfitting.
If the validation split causes its lower accuracy than 0.8 I'd look at gammas = ... 0.0001,0.001,0.1 ... because you have missed 0.01 and that will affect your model. It shouldn't give a lower score than the default parameters - it's not possible and something is wrong with the code and  application of parameters.
Once you've the initial parameters you can refine the parameter search around the initial values - nothing wrong with that and that can even be coded (automated). The only minor issue is reducing the size of the training split (60% is a bit low). That can cause issues in which case the split is 70:20:10, i.e. 10 for validation.

Thanks for your reply, but I tried what you suggested but it gave me lower accuracy values.

@user17818805 answered above below the last -------

I understood you the least lost data will infer my model. I can't improve my accuracy, I've changed the model inputs leaving only the OTUs, I've tried the PCA but to no avail.

I did'nt get that 'Once you've the initial parameters you can refine the parameter search around the initial values - nothing wrong with that and that can even be coded (automated). The only minor issue is reducing the size of the training split (60% is a bit low). That can cause issues in which case the split is 70:20:10, i.e. 10 for validation'.


text:
Scale on the y-axis in a tree rotated with 'layout_dendrogram()' using ggtree
I'm editing a time-scaled tree using ggtree. The dates are correctly adjusted when I apply the 'mrsd' command and 'theme_tree2()', but when I try to rotate the tree with 'layout_dendrogram()', the scale still remains on the x-axis with no dates. Is there a way to place the date scale on the y-axis in a tree rotated with 'layout_dendrogram()'?
# Code without layout_dendrogram()
ggtree(my_tree, mrsd='2022-12-31') %<+%
  metadato_arbol +
  geom_tippoint(size = 2) +
  theme_tree2()

# Using layout_dendrogram()
ggtree(my_tree, mrsd='2022-12-31') %<+%
  metadato_arbol +
  geom_tippoint(size = 2) +
  theme_tree2() + layout_dendrogram()

Without a minimal reproducible example this is only a guess, but try substituting theme_dendrogram() for theme_tree2() and see you still get the same output.

Maybe try the following via gridExtra,
library(gridExtra)
library(ggtree)
library(ape)
set.seed(2023-08-23)
my_tree <- rtree(30)
p <- ggtree(my_tree, mrsd='2022-12-31') %<+%
  metadato_arbol +
  geom_tippoint(size = 2) +
  theme_tree2() # + geom_text(aes(label=node))
gridExtra::grid.arrange(p, p %>% rotate(33) %>% rotate(52), ncol=2)


text:
scRNA: What are good dimensionality reduction/clustering parameters to get biologically plausible groupings?
I've got a moderately large set of PBMCs, over 1M cells. That means I can't easily do a grid search of dimensionality reduction/clustering parameters/methods. Some examples results I'm getting with scanpy default parameters include PCA -- a big, dense grouping with lobes radiating outwards, and UMAP -- "better" groupings in some sense, but many many tiny clusters of cells, more than are biologically plausible (hundreds) as well.
Basically I'm looking for best dimred/clustering practices for such a dataset, other than what's on offer from scanpy/Seurat. My metrics such as they are are that the major cell types should be clearly separated, and that the overall plot have some biological plausibility in terms of number of visually distinct clusters.
I'm aware of the Xiang et al Frontier Genetics 2021 paper, but the number of cells their example datasets use is much smaller AFAIK. They concluded that tSNE is the best choice of method according to their metrics, though UMAP also has some stability advantages. The Chari and Pachter PLOS Comp. Bio. 2023 paper raises some good points (a highly dimensional space can't really accurately have its relations preserved in 2 dimensions, and there are other issues of noise, interpretability, and so on. They recommend clustering from the higher dimensional data, which I do agree with. A different paper rec could be great.
PS: SCVI seems like a good option, I'm playing with that at the moment.

Hi M, thanks for the reply. The first paper I mentioned covers several dimensionality reduction methods and concludes that for the older data they tested, tSNE is best. The second one covers some of the theory and recommends perhaps side stepping dimred entirely.

As for metric, basically I'm looking for separation of major cell types and a plausible overall number of clusters. PCA for instance does separate some PBMC cell types out but not most. UMAP does separate out certain cell types (B cells are well separated for instance) but not others and seems to have lots of spurious tiny clusters.

Yes, I'll edit it in a bit -- will write up a less pithy short summary of the second paper.

I always wonder what 1M cells can answer what 100k cannot. Anyway, what you have to find out based on published datasets of that magnitude is which approaches scale to that number of cells. I never had to do that myself since our datasets are always way smaller. I would consider not to analyse this per-cell but rather try to find some sort of meta-cell approach where you bin cells into metas, to reduce the computational burden. I would not believe that with 1M cells one really benefits from pure single-cell analysis. Maybe some sort of clustering to get "pseudobulks" and then use these.

@M, just SCVI at the moment. @ATpoint, it's 1M+ cells across many individuals, and I do agree that cluster level comparisons are a good way to make the data tractable -- hence the question, really. Point well taken to compare to other similar datasets.

Aim It's important to define what all these calculations are doing: they're attempting to extract maximum "multi-dimensional" variance from a data set. The key thing here is that the methodology is old, but it is very well established (long before RNAseq ever existed). I get the point empirical stats is boring, but it's essential to consider whats the underbelly of the methods being used.
To answer the question directly

PCA+tSNE is worth considering
SCVI-tools is in a different league, and totalIV looks very interesting (what I'm tentatively recommending).

Note, the commercial sector is often more advanced than public sector: I developed ML methods which intersect it with classical approaches for this basic style of calculation. The scvi-tools however is "full-on" via a DL approach.
Note2, I've not read the papers cited except the scvi-tools documentation
Data set size I don't do PMBC calculations - I do perform PMBC post-hoc analytics BTW. I personally don't see the data set size as being prohibitive for an unsupervised learning calculation. Pyscan should be able to handle this.
PCA is the goto approach stage 1 - but there's a caveat, which is what you're hinting at. First, it's how much variance is trapped within each component. Scanpy will permit this calculation. What you want is >70% variance, classical stats for example PC1+PC2+PC3+PC4 > 70%. However, from past posts on this site then for PMBCs this might be unattainable, like no-can- exceed the variance threshold for lots and lots of PCs. That is an issue.
To go through your papers

Chari and Pachter PLOS Comp. Bio. 2023

Sure, I've not read the paper but on your report, simply taking PC1 and PC2 on bivariate axes is not cool at all. Very little of the variance can be represented. Nobody disputes this. Even if 10 PCs trapped > 70% of the variance - thats a lot of bivariate plots to assess, e.g. PC1 vs PC3, PC1 vs PC4 etc ... At least you've trapped a good proportion of the variance.
The problem becomes what happens when very little of the variance is captured. This is where Xiang et al Frontier Genetics 2021 are likely to be coming from. It's better at that point to use tSNE, but its not ideal.
Whats the solution?
Approach 1
In classical stats (not in RNAseq stuff) the output of PCA would be placed into tSNE and this will cause the variance spread across large numbers of PCs to be trapped by tSNE. This will give you nice clusters. I've not seen this approach at all in RNAseq - it is a common approach. The assumption is that loads and loads of PCs from the PCA trap sufficient amounts of the variance and tSNE will nicely represent that in 2 dimensions. It's possible that is what Xiang et al Frontier Genetics 2021 have done, but they may have simply gone straight to tSNE (see note2).
Approach 2
scvi-tool is impressive and I would recommend totalVI for dimensionality reduction, which they place straight into scanpy.
The method and usage (not the robustness) for totalIV scvi-tools is described here. They are deploying a "hidden Markov"/latent variable style model. See the note below on pss. They've used "CITE-seq RNA", which I don't exactly know what that is but hopefully that will remain relevant.

It's important to note scvi-tools incorporates lots of other approaches including semi-supervised learning (is that scANVI???), transfer-learning (assume its CNN-style calculations), a standard ANN type model (I think thats scVI). I'm simply proposing what I would guess would work in your case.
It is worth noting there are imperfections with the authors' methods - in general. Computer vision is a useful example, this used very similar approaches and it took years before serious levels of accuracy were achieved. Trapping absolute accuracy within a single deep learning model "off-the-cuff" I suggest is unlikely.

Semi-supervised learning its basically unsupervised learning but quantifying the weights of each feature - you can use it but I have never got good results from semi-supervised learning methods in general, this might be different though.

ps. If I've time (i.e. might not happen) I'll pop the link into the related question, where the low amounts of variance between PCs, and how to optimise that is partially answered.
pss. Overall scvi-kit for the models outlined need published demonstration they are robust. This part of their work appears partially missing, basically there appears a lot more functionality than published papers. I assume their other stuff is working its way through peer-review. The totalIV approach: I performed a brief search but was unable to find anything published - whilst I'm recommending it - it appears untested. What I believe is they've stolen a march on the rest of the field in general, so I'm making a best guess.


text:
Snakemake: Re-read the samples.tsv file midway through the workflow
All, I have a workflow that in some rules requires the sample's read-lengths. For this I always ran a script that would parse all fastq files and would add the read lengths to the sample table. Rules could then take the read lengths for each sample from there during execution.
However, I recently made a workflow/pipeline that downloads files from SRA, so I need to compute the read lengths at a later stage during workflow execution (after prefetch/fasterq-dump).
I added a step to update the sample table, however, I think the sample table is only read at the very beginning of the workflow.
What would be the nicest way, the Snakemake preferred way, to re-read the samples.tsv file before a rule is executed?
In general I like storing meta data in the sample table so, I'll probably use this mechanism more often in the future, if I find a nice methods to re-parse the samples table.
Right now I do this in rules:
rule deletion_profile:
    input:
        rseqc_input_bam,
        rules.add_max_read_length.output.done_file
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    params:
        read_length = lambda wildcards: samples.loc[wildcards.sample, 'read_length']
    shell:
        '''
        deletion_profile.py -i {input} -l {params.read_length} -o {rseqc_dir}/{wildcards.sample}

Thanx.

You're right, it doesn't re-run things again for each job, so your original definition of samples is still the one it's working with.
It seems like it goes against the grain of Snakemake to overwrite a file in-place that acts as input to a rule, so I think the most Snakemake-y way you could do it would be to treat your samples.tsv as part of the input/output graph like everything else.  You could have a rule with samples.tsv as input and samples.somethingelse.tsv as output where that latter one has read lengths and whatever else gets automatically added.  (Alternatively, you could just have your updated samples.tsv read within the rule itself rather than relying on a shared copy loaded earlier, which is kind of halfway to the above idea.)
Either way, maybe something like:
rule deletion_profile:
    input:
        bam=rseqc_input_bam,
        samples="samples.extended.tsv"
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    run:
        samples = whatever_loads_samples(input.samples)
        read_length = samples.loc[wildcards.sample, 'read_length']
        shell("deletion_profile.py -i {input.bam} -l {read_length} -o {output}")

(On a side note, you can use {output} rather than repeating yourself for that part.)


text:
Can not launch bcftools using python's subprocess module, as it only accepts first command of commands list
I am trying to remove samples from a chromosome vcf file. I wrote a function that takes chromosome number and a list of samples to remove. When I try to run bcftools using subprocess module it only runs bcftools, as if I was running subprocess.run([bcftools], shell=True, env=my_env) instead of:
subprocess.run(['bcftools', 'view', in_vcf, '--samples', drop_samples, '-Oz', '>', out_vcf], shell=True, env=my_env)

The function:
def remove_individuals(chr_samples):
    my_env = os.environ.copy()
    my_env['PATH'] = f'/home/aydar/bin/bcftools_1.11/bin/:{my_env["PATH"]}'
    def remove_indiv(chr, indiv):
        in_vcf = f'data.filtered-updated-{chr}.vcf.gz'
        out_vcf = f'data.filtered-updated-{chr}-faulty-samp-rem.vcf.gz'
        drop_samples = "\'^"+",".join(indiv)+"\'"
        proc = subprocess.run(
            ['bcftools', 'view', in_vcf,
             '--samples', drop_samples, '-Oz', '>', out_vcf],
            shell=True,
            env=my_env
        )
        print(proc.args)
    chr_samples.apply(lambda x: remove_indiv(x.index[0], list(x.values[0])))

function input:
13,[2500000075_2500000075,2000000171_2000000171,1800000179_1800000179,330134048_330134048]

When I run the same command in the terminal everything works.
The command I test with in terminal:
bcftools view data.filtered-updated-13.vcf.gz --samples '^2500000075_2500000075,2000000171_2000000171,1800000179_1800000179,330134048_330134048' -Oz > out_test_vcf.gz

Thanks for your question and answer. Using popen is the standard approach here than using run. Anyway could you very kindly mark your answer as "accepted"? This help site stats and prevents the question from being re-introduced as unresolved later.

Never mind, I figured it out. As it turned out I had to write to file like so and extra quotes ("\'") around samples were not needed. Also shell=True does not work well with lists, only strings:
def remove_individuals(chr_samples):
    my_env = os.environ.copy()
    my_env['PATH'] = f'/home/aydar/bin/bcftools_1.11/bin/:{my_env["PATH"]}'
    def remove_indiv(chr, indiv):
        in_vcf = f'data.filtered-updated-{chr}.vcf.gz'
        out_vcf = f'data.filtered-updated-{chr}-faulty-samp-rem.vcf.gz'
        drop_samples = "^"+",".join(indiv)
        out = subprocess.run(
            ['bcftools', 'view', in_vcf,
             '--samples', drop_samples, '-Oz'],
            env=my_env,
            capture_output=True
        )
        with open(out_vcf, 'wb') as f:
            f.write(out.stdout)
    chr_samples.apply(lambda x: remove_indiv(x.index[0], list(x.values[0])))

Remember to use 'wb' if you like me and use Oz output option, because it returns bgzipped file.

thanks for accepting your answer and please feel free to ask a question any time here.


text:
Variation in 16S rRNA between assemblers - how do I know which is real?
I have a low-diversity metagenome (~11 bins > 80% completion). Out of the bins, 3 are of interest to me. None of the lower-completion bins that can be identified are from the group of interest. So let's call the bins of interest 1, 2 and 3.
The raw reads had been assembled separately with MEGAHIT and SPAdes (meta). Out of them, only the SPAdes assembly was binned / refined. There are 3 16S rRNA in the MEGAHIT assembly (x, y, z) and 3 16S rRNA in the SPAdes assembly (a, b, c).
Bin 1 is associated with sequence 'a', which is identical to 'x'.
Bin 2 is linked with SPAdes sequence 'b', which has no MEGAHIT equivalent. Bin 3 is not associated with any 16S rRNA (even with MarkerMAG).
Sequences 'c' and 'z' are very similar. Sequence 'y' is an outlier.
There are small hints from phylogenetics that bin 3 could be associated with sequences c/z. But then, could sequence 'y' just be a misassembly? Or is it that somehow my final bins are missing an entire genome of interest? If the 16S rRNA for 'y' is real, I would expect a genome divergent enough from the others that it wouldn't get binned together with them. This is what I am trying to understand. Would you happen to have any advice / suggestions on what I could to to test this?
Initial approach: assessing scaffold quality with the BBMap stats.sh . It's a bit difficult to tell, since SPAdes has about 10x more contigs. Worse N50/N90 but more long contigs and a better 'general assembly score'.
Additional thoughts: The new bins seem too divergent from known genomes to use a reference for misassembly estimation.
Thank you for your time!

This is easy @Laura - I'm a bit surprised by the results Spades meta didn't have a good reputation, but the criteria assessed was not what you are looking for.
Anyway take all outputs and place them in large multi-species alignment,

e.g. for each output perform a Blast for maximum genetic diversity for each sequence and sample across the large Blast  output. Take note if there are multiple partial hits against diverse bacteria (online Blast has a diagrammatical summary).
You have a good understanding of your taxonomic group and intersperse that with other stuff, i.e. fill in the diversity.
Assemble it, e.g. muscle5 superfast option.
Perform a split decomposition phylogeny or equivalent which identifies phylogenetic networks.

SplitsTree is the historic package for this and it turns out that split decomposition was not the best, I think "neighornet" was superior. What this is does is identify mosaic gene signal in your data. It provides a visual representation. Thus if your data is a artifactual mix of two or more taxa this will cause a network effect between the two parents and your newly mis-assembled sequence will sit between the two in a two "diamond same". This is the antonym of a bifurcating phylogeny.
The final step is the assembly that produces sequences that minimises phylogenetic networks is the true sequence. If you've not encountered split-decomposition/networks before a verbal explanation might seem confusing but when you run through a few examples and try it on your data you'll get the grasp.

Much appreciated! Good idea. Tested. Couldn't see obvious signs of a chimera, but the alignment showed that most of the short 'outlier' sequence's distinguishing features were not present in any other (even distantly related), including a lot of weird Gs. Chimeric with something else maybe? Curious thing is, the depth of that MEGAHIT contig seemed to be pretty good (~120x). I had thought I'd done my homework by quality-controlling the reads with BBDuk, but it's a bit unsettling to see how software can be so certain and so wrong at the same time. Should look more into missassembly detection...

It tentatively sounds like you might have found a new species. Just keep in mind the diversity of bacteria in the analysis affects the ability to detect a possible chimera, i.e. if the "parents" are not present it could throw a false negative -------- Weird "Gs" high AG are a known feature of bacteria: particularly if its genetically distinct. Anyway if Splitstree if done correctly it's a sensitive method for chimeras: but you might have to shift the sample size - remove distant sequences which can cause a collapse in networks.  Megahit is respected BTW.

It's more like...The positions in questions are 100% conserved within related bacterial orders, and here's this short sequence with ~9 of these otherwise conserved positions being different. I guess I can't prove it's not real (part of my issue...), but it doesn't look...correct?

What you are saying sounds perfectly reasonable - keep in mind I can't see the data nor analysis. Proof ... it's the first stage requiring downstream analysis, so either leave it as a "putative result" or go into the complexity of metagenome assembly by assembling other single copy gene targets like atpD.


text:
Determining fragment mean and fragment stdev for MaSuRCA config file
Similar to this unanswered question on Biostars, I am using MaSuRCA for the first time and want to know how other MaSuRCA users are determining fragment mean and fragment stdev. My understanding is that FastQC measures read length, which differs from fragment length.
This paper that used MaSuRCA for a bird genome assembly mentions the parameters they used but not how those parameters were determined: "MaSuRCA was run applying the following parameters: fragment mean (422), fragment stdev (312), and estimated genome size (1.2 Gbp)."
I am using Illumina, Nanopore, and PacBio reads for fungal genome assembly. Any advice on which tool(s) to use to determine fragment mean and stdev for MaSuRCA would be appreciated.

@M__ Added the quotation and clarified that I am doing fungal genome assembly

Thanks for the edits and added information @juliadouglasf .. difficult genome BTW (Mycorrhizae I assume), but definitely an interesting and worth while area and good luck with the question.


text:
Gene Set Comparison Without Expression Data
I have been looking all over the web to find some answers to my problem but unfortunately, I was unsuccessful. I wish to determine whether an a priori defined set of genes in my case genes associated with Epidermolysis bullosa shows statistically significant, concordant differences between two biological states (e.g. phenotypes).
I have a list of genes from my own experiment and a predefined list of genes associated with Epidermolysis bullosa, and want to find out if the genes of Epidermolysis bullosa are overrepresented in my list.
Here is an article about GSEA, which is related to what I want to do:
https://www.pnas.org/doi/epdf/10.1073/pnas.0506580102

so I have the following data:
I'm researching patients with severe or mild disease; in this research, I have found a list of  200 genes after conducting a SKAT test in R (https://www.rdocumentation.org/packages/SKAT/versions/2.2.4/topics/SKAT) with a p-value <0.05 ( which I presume can cause a severe disease instead of mild  ) I want to focus on this 200 significant genes that can perhaps affect the disease level, this genes can be  -> the L list as described in the article above.
Also, the S list of predefined genes ( as is called in the article above )is associated with Epidermolysis bullosa disease. and I want to know if the genes from the Epidermolysis bullosa are overrepresented in my L list.
After reading the article and many more resources, I understood that to perform GSEA  (https://www.gsea-msigdb.org/gsea/index.jsp), I had to have gene expression data which  i don't have .
so I wish to understand if there is a tool that could do a gene set comparison without having the expression data whether it's in R or Python that uses my gene list and predefined list ( in my case associated with Epidermolysis bullosa)  ,thank you :)

You can't do a GSEA on a set alone because the differentiating statistic (or, occasionally, the rank) is used to calculate the enrichment score that is produced by the GSEA algorithm.
The other gene list type tools that I'm aware of (that my biologist colleagues use) are DAVID and PANTHER. But those are general discovery algorithms, not for comparing specific lists.

thank you , is there a way to rank the genes not based on expression data ? Or a tool that calculates the enrichment not using the   GSEA  algorithm?

@ElizaR Sounds a bit like a pathway analysis (e.g. see gprofiler). Perhaps a hypergeometric test


text:
GenomicsDB and recalibration recommended for cohort of mixed exome (GATK)
We have a cohort of around 3k exomes, most of them from different capture kits, such as Agilent v6, v8, Nextera, Twist, Clinical Research exome and many others. Right now we are creating genome db on a per chromosome basis and combining the vcf.
We were wondering the following questions -

Will the mix of capture kits introduce more false positive variant calls in the recalibration step

If we restrict the genomics DB to only include regions covered by capture kit, will we see a performance boost (Ex. instead of entire chr2, we restrict it to regions covered by all capture kits).


text:
Why does my HWE QQ plot have extreme deviation and what does it mean?
This question was also asked on Reddit
I have recently completed my thesis and one of the comments was that I report on why this graph looks this way. I have tried to find a reason but the closest I can come to is that it means that there are false positives? Or that this is possibly due to population stratification?
The trait of interest was alcohol dependency with covariates of sex and age, in a non african american population.
I am also unsure as to how this will impact my research or the implications of this?

I don't understand why this population (non-African American) "reacted" this way - do I provide a different null model? The previous populations had "relatively normal" plots. They were all run through the same pipeline.
How the data was calculated:

There were 1061 individuals in the non-African American dataset. At QC phase 0, there were 2 individuals who had discordant sex information and 87 duplicate SNPs. When the MAF value was set to 0.000015, 125708 SNPs were removed when the MAF. Both the mind and geno parameters were set to 0.2, resulting in 127357 and 0 SNPs to be removed respectively. When setting the HWE to $1 \times 10^{-12}$, 4495 SNPs were removed. After all these parameters were set and the pipeline completed, there were 851773 SNPs left from 1044 individuals.

This much deviation from the expected distribution looks like a modeling issue, i.e. that the statistical model use to represent the null hypothesis is incorrect.
I can't see any values that lie along the expected line. It's more common to see a generally linear trend, with a few outliers just at the end.
It could very likely be due to population stratification. You have defined this population as "non-African American", which seems to me to be quite a vague definition that would capture a lot of different population groups with distinct genetic histories.
In my own research many years ago, I found substantially different haplotypes within the ADH gene cluster between Māori and European individuals, and I could imagine that not adjusting for population structure in other populations would have an impact on association scores.

The level of discordance is extreme, particularly for humans. This is either a publication in Nature/Science and massive advance in the field (especially humans), or the genetic factors have massive confounders.
In theory, I'm aware there must be a genetic component to alcoholism, there is a biochemical component about the rate of processing alcohol toxins. However, non-genetic (especially social/community) affects of alcoholism could predominate.
What I personally think is happening is there are very strong non-genetic confounders which are disrupting the model and this is getting confused with neutral HWE. For example, if alcoholism predominates in a community or social group the transmission is social not genetic, so the confounder is extreme because they are genetically related, but the gene(s) are not being inherited. In phylogenetic terms it will make the "gene" appear far higher up (far deeper) in the tree than it actually is. The "lone" drinker is probably the model, but group behaviour and peer-pressure is probably more realistic. Equally community effects could be suppressing alcoholism, such as religious belief.
Thus if a non-genetic model of alcoholism, (loads of social factors) resulted in a stronger fit by for example 3 logs difference, whilst the genetic model is 50 logs difference, then that is the answer and points to the need for careful population sampling for GWAS. Upvotes J Human Gen is cool.

As a disclaimer, I do not understand not human genetics outside immunology. I do understand multiple linear regression models and failing to account for a pivotal factor(s) causes unusual effects. An obvious example in this case is religious belief: thus alcoholism genes could be present but the environment is not presented. For example, you might conclude long hair caused a reduction in human growth because they tightly correlated in humans. However, once gender is a factored into the model, the result is clear hair-length associates with gender not height.


text:
How can I retrieve the GRCh38 coordinates of a list of rsids?
I have a list of about 100,000 rsids and I want to get their genomic coordinates on the GRCh38 genome build. Is there a command line tool that allows me to do this? If yes, which one?
I have tried bcftools but, given the error message I got, I believe that it does not support search by rsid:
bcftools view -R rsids_chr1_subset.txt -f PASS -v snps dbSNP_v156_GRCh38.p14.vcf.gz > rsids_chr1.vcf
[E::bcf_sr_regions_init] Could not parse 1-th line of file rsids_chr1_subset.txt, using the columns 1,2[,-1]
Failed to read the regions: rsids_chr1_subset.txt

The tricky thing is doing this fast. If I remember correctly, I have successfully used an implementation written in C a long time ago that did just this. It also had a Python version. The name escapes me at the moment. However, this just reminded me that I ended up writing my own implementation because there was an issue. I will make an answer.

Several years ago I created rsidx to address this particular use case. https://github.com/bioforensics/rsidx
# VCF should be sorted by genomic coordinates and indexed by tabix
rsidx index dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx
rsidx search dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx rs3114908 rs10756819

rsids_chr1.txt contains a list of rsids, one per line.
bcftools view --include ID==@rsids_chr1.txt dbSNP_v156_GRCh38.p14.vcf.gz -Ou | bcftools query -f '%CHROM  %ID %POS  %REF  %ALT{0}\n' > dbSNP_v156_GRCh38.p14_AF.vcf

This was part of a custom liftover function which handled multiple alleles and reverse compliments. I do not know if it will be useful. I believe there was a separate thing somewhere with lookup that I may add. It handled apparent non-1:1 mappings in addition to platform-specific "rsIDs."
converter = get_lifter(current, target)
def convert_coordinates(input_file, output_file, converter):
        line_count = 0
        skipped = 0
        strand_flips = 0
        with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
            for line in infile:
                # Skip headers
                if line.startswith("#"):
                    outfile.write(line)
                    continue

                # VCF files are tab-separated
                parts = line.strip().split('\t')

                # Get the chromosome and position
                chrom = parts[0]
                pos = int(parts[1])

                # Convert the coordinates
                try:
                    new_coords = converter.query(chrom, pos)
                except KeyError:
                    #print(f"Chromosome {chrom} not found in the liftover chain file.")
                    skipped += 1
                    continue

                # If there's no new coordinates, we skip this variant
                if not new_coords:
                    skipped += 1
                    continue

                # Get the new chromosome, position and strand
                new_chrom, new_pos, strand = new_coords[0]

                # If the new strand is negative, reverse complement the REF and ALT alleles
                if strand == '-':
                    parts[3] = reverse_complement(parts[3])
                    parts[4] = reverse_complement(parts[4])
                    strand_flips += 1

                # Replace the old chromosome and position with the new ones
                parts[0] = new_chrom
                parts[1] = str(new_pos)

                # Write the line back into the file
                outfile.write('\t'.join(parts))
                outfile.write('\n')

                # Print a progress statement every 1000 lines
                line_count += 1
                if line_count % 100000 == 0:
                    print(f"Processed {line_count} lines")

        # Calculate and print the percentage of variants skipped
        percentage_skipped = (skipped / line_count) * 100
        percentage_flips = (strand_flips / line_count) * 100
        print(f"\nFinished processing. Out of {line_count} total lines:")
        print(f"{percentage_skipped}% of variants were skipped due to missing chromosomal data.")
        print(f"{percentage_flips}% of variants required a strand flip.")


text:
Help with error DiffBind
This question was also asked on Biostars and the Bioconductor Forum
I want to find differential accessibility between conditions (WS, WF, MS, MF). I can compare any pair conditions without error except this pair MS vs MF. Data is the bam and broadpeak files from ATAC seq pipeline and the code is like a pipeline of the DiffBind R package. The sample is human and uses Illumina sequencing technology.
Would you please suggest what wrong in this case? I pretty much run the same code as I used for another data set and didn't have this error. Thank you so much!
result <- dba(sampleSheet=samples)
result <- dba.blacklist(result)
result <- dba.count(result, bParallel=FALSE)
result <- dba.count(result)
result <- dba.normalize(result)
result <- dba.contrast(result, minMembers=2)
result <- dba.analyze(result)
result <- dba.plotProfile(result)

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

There is post about this error but I still got error when trying those solutions. Such as the order of the function, the path to the bam and broad peak files in sample sheet.
dba.show(result, bContrasts=TRUE)
     Factor    Group Samples  Group2 Samples2 DB.DESeq2
1 Condition diseased       2 control        2         0

I don't know why the value of DB.DESeq2 in this case was 0 which maybe the reason of the error.
This is the result after running dba.contrast:
4 Samples, 89585 sites in matrix:
   ID Tissue   Factor Condition Replicate    Reads FRiP
1  MS    non ATAC-seq   control         1 23731398 0.17
2 MSR    non ATAC-seq   control         2 10499000 0.16
3  MF    non ATAC-seq  diseased         1 31535510 0.18
4 MFR    non ATAC-seq  diseased         2 18488559 0.18

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

If I try to compare MS vs MF, I will get that error but if I remove MS and replace by other sample, I don't have that error anymore.
I don't have that error if I add more samples and the value of DB.SESeq2 at the first row from dba.show() is not 0.
Full samplesheet:
SampleID,Tissue,Factor,Condition,Treatment,Replicate,bamReads,ControlID,bamControl,Peaks,PeakCallers
ws1,non,ATAC-seq,ws,,1,ATAC_seq_output/bwa/merged_library/ws_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP1.mLb.clN_peaks.broadPeak,bed
ws2,non,ATAC-seq,ws,,2,ATAC_seq_output/bwa/merged_library/ws_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP2.mLb.clN_peaks.broadPeak,bed
wf1,non,ATAC-seq,wf,,1,ATAC_seq_output/bwa/merged_library/wf_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/wf_REP1.mLb.clN_peaks.broadPeak,bed
wf2,non,ATAC-seq,wf,,2,ATAC_seq_output/bwa/merged_library/wf_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/wf_REP2.mLb.clN_peaks.broadPeak,bed
ms1,non,ATAC-seq,ms,,1,ATAC_seq_output/bwa/merged_library/ms_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ms_REP1.mLb.clN_peaks.broadPeak,bed
ms2,non,ATAC-seq,ms,,2,ATAC_seq_output/bwa/merged_library/ms_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ms_REP2.mLb.clN_peaks.broadPeak,bed
mf1,non,ATAC-seq,mf,,1,ATAC_seq_output/bwa/merged_library/ws_REP1.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP1.mLb.clN_peaks.broadPeak,bed
mf2,non,ATAC-seq,mf,,2,ATAC_seq_output/bwa/merged_library/ws_REP2.mLb.clN.sorted.bam,,,ATAC_seq_output/bwa/merged_library/macs2/broad_peak/ws_REP2.mLb.clN_peaks.broadPeak,bed

Hi @gringer, I added the sample sheet in the question. Any other pair comparisons have non zero DB.SESeq2 value.

The conditions in that sample sheet (ws / wf / ms / mf) don't match the dba.contrast result (control / diseased). Can you show the result of dba.contrast after doing another pair comparison?

It is the full sample sheet. Because I don't want to compare ws and wf so I delete them. Maybe there is a way to manipulate the output but I still haven't known so I delete ws and wf.

How did "diseased" and "control" get into the first result after dba.contrast? They don't appear on the sample sheet at all. Why does the ID column look different from the sample sheet? I'm being pedantic about this because it seems important for the specific error that you're having.

I replaced ms in condition column by control and mf by diseased. Thank you. Hope it is clearer now.

I have very little idea what's going on here, but in the absence of any other answers or help, I'll try to work through this from a no-knowledge start point to see if I can hopefully help shift this problem closer towards a solution.
The specific error you have encountered is No valid contrasts/methods specified. Looking at the DiffBind Vignette, it provides this information about defining contrasts:

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

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

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

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

Specifically relating to MS vs MF, I can see why there would be an issue with the condition contrast, especially if it's using DESeq2 in the background. The short answer is that there is not enough variation within the groups.
DESeq2 cares about making sure that conditions added to the model are separable from samples. This is further explained in the DESeq2 manual here. The issue is that there's no way to distinguish between the control condition and the MS sample (or, similarly, between the diseased condition and the MF sample). In order to get past this error, you need both control and diseased in your MS group or in your MF group (or both, if you want to be really rigorous), and ideally a few more replicates of each condition.


text:
How to rearrange gene ontology obo file?
I have a gene ontology obo file, and would like to rearrange the obo format file in a single column to a text file in three columns comprised of 'parent (id:GO:xxxxxxx)', 'child (is_a: GO:xxxxxxx, part_of GO:xxxxxxx, or positively/negatively_regulates GO:xxxxxxx. For example...
Thank you for your help in advance!
Original file:
[Term]
id: GO:0000001 # parent
name: mitochondrion inheritance
namespace: biological_process
def: "The distribution of mitochondria, including the mitochondrial genome, into daughter cells after mitosis or meiosis, mediated by interactions between mitochondria and the cytoskeleton." [GOC:mcc, PMID:10873824, PMID:11389764]
synonym: "mitochondrial inheritance" EXACT []
is_a: GO:0048308 ! organelle inheritance # child
is_a: GO:0048311 ! mitochondrion distribution # child

[Term]
id: GO:0000015 # parent
name: phosphopyruvate hydratase complex
namespace: cellular_component
def: "A multimeric enzyme complex, usually a dimer or an octamer, that catalyzes the conversion of 2-phospho-D-glycerate to phosphoenolpyruvate and water." [GOC:jl, ISBN:0198506732]
subset: goslim_metagenomics
synonym: "enolase complex" EXACT []
is_a: GO:1902494 ! catalytic complex # child
relationship: part_of GO:0005829 ! cytosol # child

[Term]
id: GO:0000122 # parent
name: negative regulation of transcription by RNA polymerase II
namespace: biological_process
alt_id: GO:0010553
alt_id: GO:0045816
def: "Any process that stops, prevents, or reduces the frequency, rate or extent of transcription mediated by RNA polymerase II." [GOC:go_curators, GOC:txnOH]
synonym: "down regulation of global transcription from RNA polymerase II promoter" RELATED []
synonym: "down regulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "down-regulation of global transcription from RNA polymerase II promoter" RELATED []
synonym: "down-regulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "downregulation of global transcription from RNA polymerase II promoter" RELATED []
synonym: "downregulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "inhibition of global transcription from RNA polymerase II promoter" RELATED []
synonym: "inhibition of transcription from RNA polymerase II promoter" EXACT []
synonym: "negative regulation of gene-specific transcription from RNA polymerase II promoter" RELATED []
synonym: "negative regulation of global transcription from Pol II promoter" RELATED []
synonym: "negative regulation of transcription from Pol II promoter" EXACT []
synonym: "negative regulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "negative regulation of transcription from RNA polymerase II promoter, global" RELATED []
is_a: GO:0006357 ! regulation of transcription by RNA polymerase II # child
is_a: GO:0045892 ! negative regulation of DNA-templated transcription # child
relationship: negatively_regulates GO:0006366 ! transcription by RNA polymerase II # child

Desired format:
Parent       Child       Relationship
GO:0000001   GO:0048308  is_a
GO:0000001   GO:0048311  is_a
GO:0000015   GO:1902494  is_a
GO:0000015   GO:0005829  part_of
GO:0000122   GO:0006357  is_a
GO:0000122   GO:0045892  is_a
GO:0000122   GO:0006366  negatively_regulates

Thank you for the advice. I've just reformatted the example file.

One potential solution using awk:
awk 'BEGIN {
    OFS = "\t"
    print "Parent", "Child", "Relationship"
}

/GO:[[:digit:]]{1,}/ {
    if ($1 ~ /^id:/) {
    parent_id = $2
} else if ($1 ~ /^is_a:/) {
        gsub(":", "", $1)
    print parent_id, $2, $1
} else if ($1 ~ /^relationship:/) {
        print parent_id, $3, $2
    }
}' test.txt
Parent      Child       Relationship
GO:0000001  GO:0048308  is_a
GO:0000001  GO:0048311  is_a
GO:0000015  GO:1902494  is_a
GO:0000015  GO:0005829  part_of
GO:0000122  GO:0006357  is_a
GO:0000122  GO:0045892  is_a
GO:0000122  GO:0006366  negatively_regulates

Including comments:
awk 'BEGIN {                      # before processing the first line
    OFS = "\t"                    # set the output field separator to tab
    print "Parent", "Child", "Relationship"    # and print the header
}

/GO:[[:digit:]]{1,}/ {    # for any line that contains a GO term
    if ($1 ~ /^id:/) {            # if the line starts with id:
    parent_id = $1                # set the parent id
} else if ($1 ~ /^is_a/) {        # then, for lines that start with is_a
        gsub(":", "", $1)         # remove the colon
    print parent_id, $2, $1       # and print the parent, child and relationship
} else if ($1 ~ /^relationship:/) {  # then, for lines that start with relationship
        print parent_id, $3, $2   # print the parent, child, and relationship
    }
}' test.txt

You're welcome - I'm glad you solved your problem @ojikj2023 - but comments are not the right place to say 'thank you'; please see What should I do when someone answers my question?.

It's my first time to post a question and get help. I've just done voting following your comment!


text:
How do I delete an amino acid residue from a loop and then join the loose ends?
I want to delete an amino acid residue from a loop and then join the loose ends. I intend to increase the thermal stability of the protein by reducing the size of the loop and its agitation. Pymol just deletes the residue and breaks the loop.
I intend to increase the optimal temperature of the protein, even if it costs some reduction in activity. The objective is to reduce protein agitation due to temperature increase. I will try point mutations that increase hydrogen interactions and eventually delete a loop located glycine to reduce its length and consequently its thermal agitation. All this work is aimed at preserving the catalytic site, but some loss of activity is acceptable.
As a complement to my question, I present an excerpt from the paper by Mahdie Rahban et al (2022) ; Thermal stability enhancement: Fundamental concepts of protein engineering strategies to manipulate the flexible structure.

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

Does anyone have a suggestion how to solve this problem? I already used Modeller to insert an amino acid residue in place of a deleterious mutation, could Modeller be used to reconstruct the loop? How do I do it?

Modeller .. no thats for point mutations

https://salilab.org/modeller/9.15/manual.pdf 2.3 Loop optimization pag 26

It depends on what software you have in hand and want to use. But basically, when you delete/add a residue, you can't just stick two parts together since it affects the whole conformation.
I suggest using AlphaFold2, and input it with different sequences (including deletion mutation) to see how it affects the structure.
Here is a link to AF2 colab:
https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb


text:
Getting rid of duplicates in a dictionary
I have a tsv file that lists the reads and read lengths from a FASTA file but some reads are duplicated - that's just from the analysis I did previously - but I want to only take one instance of the read and have it be the SHORTER length than the other duplicate.
For instance, here is the beginning of the tsv file:
read                                    length
91290b97-5ca1-489d-b488-17afc1d217e7    3584
dad08436-3cd6-40a1-9af1-b846bd24a815    3787
f4a7bbdf-de62-48a8-a8b0-06076e251820    3637
f4a7bbdf-de62-48a8-a8b0-06076e251820    4717
21ddb2f2-654e-4124-8226-91c5a8d11a97    3786
6a2de460-a661-4c9e-ac1d-adfc3e5841e2    3520

You can see that f4a7bbdf-de62-48a8-a8b0-06076e251820 is found twice with different lengths 3637 and 4717. I want to go through and take the shorter length of that read and include it in the final tsv file, resulting in this structure:
read                                    length
91290b97-5ca1-489d-b488-17afc1d217e7    3584
dad08436-3cd6-40a1-9af1-b846bd24a815    3787
f4a7bbdf-de62-48a8-a8b0-06076e251820    3637
21ddb2f2-654e-4124-8226-91c5a8d11a97    3786
6a2de460-a661-4c9e-ac1d-adfc3e5841e2    3520

I'm trying to do this in python and have started with making the tsv file a list:
newlist = []
duplist = []

with open(tsv_path) as file:
    tsv_file = csv.reader(file, delimiter="\t")

    for line in tsv_file:
        if line[0] not in newlist:
            newlist.append(line[0])
        else:
            duplist.append(line)

    print("All the duplicate reads: ", duplist)
    print("The normal reads: ", newlist)

Which works to contain the first instance of a read in one list and the second/third instance of the read to go in a separate list with their lengths. However, the shorter length of a read is not always the first instance of a read. Do I need to include another loop that looks at the lengths and compares each instance of a read? Is there an easier way to do this?

Not sure how you'd do this in python but you need a group-by operation here. Group by the read name and pick the min length per group.
EDIT: With pandas, you can do:
df.groupby(['read'])['length'].agg('min')

That's so much simpler and it worked! Thank you!

You can do this with the sort command-line tool. This command pipe sequence sorts numerically by the second field, then chooses unique elements based on the first field
sort -r -k 2n,2 read_lengths.txt | sort -r -u -k 1,1 > unique_lengths.txt

[the additional -r bits are not strictly necessary; they help make sure that the headings stay at the top]


text:
nextflow: tumor normal sample - how to make code organic
I'm developing pipeline in nextflow it uses tools: fastp, bwa mem index, bwamem, gatk mark duplicates, gatk setupnmd, gatk applybqsr, gatk recalibrate. I do this using FASTQ.gz of normal and tumor DNA samples.
I'm interested to use mutect2, lancet and strelka in nextsteps. However, I am unable to write code in the pipeline that makes sense. I use loads of shell scripting in mutect2 to use tumor-normal.
To expand, I've input in FASTQ as tiny_n, tiny_t, where _n is normal and _t is tumor. I run pipeline until applybqsr for each individually. However, as I reach mutect2, I check if sample contains _n or otherwise, I then create an empty file for the output as required in the mutect2 process.
How do I overcome creating dummy/empty file?
Please see below code of main.nf:
params.raw_files="/sc/arion/projects/tiny/all_tumors_normals/tiny_*{R1,R2}.fastq"

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

        read_pairs_ch = Channel.fromFilePairs(params.raw_files)

        FASTP(read_pairs_ch)

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

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

mutect2(applybqsr.out.recal_bam)
}

below is code of mutect2 to handle only one sample, run for normal, create dummy/empty file for tumor sample, mutect2_process.nf:
process mutect2 {
        maxForks 3
        debug true
        tag "${sample_id}"
 memory '15 GB'

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

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

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

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

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

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

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

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

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

        fi

        """
        }

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

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

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

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

    ...

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

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

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

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

        }
        .set { cases_ch }

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

    mutect2( joined_ch, ref_fasta, targets )
}

process mutect2 {

    tag { sample_id }

    module 'gatk/4.1.3.0'

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

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

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

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

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

I really don't understand what you did with .branch in applybqsr. Nextflow is so mysterious. How are you defining normal: and tumor: without define keyword? Also, I have samples as: MM_tiny_n_R1.fastq.gz MM_tiny_t_R1.fastq.gz Thank you for your reply.

Also, where is cases.tumor stored/created?

@DeathMetal Sorry that should have been cases_ch.tumor. Fixed now. The boolean expressions we define in our branch closure can be given unique labels, which we can then use as named channels.

@DeathMetal With your current glob pattern, the fromFilePairs factory method should give you the sample names: MM_tiny_n and MM_tiny_t. I assumed your FASTP and dowstream processes also return tuples containing the sample names/ids, but I might be wrong..

This is complicated for me, sorry. Would you mind explaining  in tumor: sample_id.endsWith('_T')  return tuple( case_id, sample_id, bam ) what is case_id? Second, how does tuple has three elements in it?

If you get a chance, I've a new question at https://bioinformatics.stackexchange.com/questions/21310/nextflow-missing-value-declared-as-output-parameter-output-variable-defined

@DeathMetal The case_id is really just the patient or case or donor from which the samples are derived. For example, MM_tiny would be the case or donor id here. Tuples can contain any number of elements. The elements don't even need to be of the same type necessarily.


text:
Compare in Python / PyRanges two genomic intervals to determine whether they're separate, partially overlapping, or contained
As part of a larger analysis, I need to evaluate pairs of genomic regions (regionsA and regionsB) and perform different calculations for each pair depending on whether:

They're completely separate (no overlap)
regionsB partially overlaps regionsA (start or end of regionsB does not overlap; both are handled separately)
regionsB is fully contained in regionsA

At this time, regionsA and regionsB are represented as Pandas DataFrames. Likewise, I'm doing all possible combinations:
if (region_b.start < region_a.start and
    region_b.end < region_b.start):
[...]
elif (region_b.start < region_a.start and
     (region_b.end > region_a.start
     and region_b.end < region_a.end )):

Given that there are many conditions I check for, this makes the code long, a little hard to read, and too verbose. Is there anything I can use (without making the code harder to read) to make these checks simpler, something like PyRanges (or PyRanges itself, although I don't think there's API for it)?
If not, I'll just live with what I have.


text:
Finding linkage disequilibrium when I have a kinship matrix but not pedigree
I have a bunch of samples of a non-model organism divided in two groups: parental generation and offspring. I do not have a pedigree so I do not actually know who is a parent of who, and which are the sibships. I have a bunch of SNPs that I can use to create a kinship matrix, but I cannot resolve families just using it.
I want to perform linkage disequilibrium (LD) between the SNPs. I can calculate the correlation between them, but I feel I am losing a lot of power because there is kinship between samples. Is there a program that could take the kinship matrix in account when calculating the LD?
I could, in theory, reconstruct the pedigree from the markers. I run multiple instance of Colony but the results are underwhelming -- most samples are not allocated (admittedly this raises questions on the SNP situation, but this is a whole different issue).

Welcome to the site. Why not just do LD directly?

what do you mean?  I have only a vague idea of actual distances (and surely no genetic distance) because it is a non model organism. All I can do is a correlation

Well given my data -- unphased, no pedigree, it's not I can get much more than a correlation, hence my question -- I can create a kinship matrix, and maybe use it a prior to have some idea of the segregation and the LD, but I have no idea if anyone made some software able to do that

yes, LD is linkage disequilibrium

I originally thought you were working on humans. I spotted "non model organism", that changes things alot, albeit I guess humans are a non-model organism in a way. The question is what the organism. LD has been applied to every kingdom of life, except viruses. Arlequin isn't bad for relatively small to medium data sets (but tricky to get the format right - thats a separate question). If it's a genome sequence then I still suspect its plink. There are a number of recent LD programs BTW.
I did notice Colony was written in Fortran! Wow, the author must be a physicist at ZSL anyway...
If you have a defined number of SNPs - sounds small to medium data set - then I'd use Arlquein http://cmpg.unibe.ch/software/arlequin35/Arl35Methods.html, Arlquin does loads of pop gen BTW its very famous for AMOVA BTW (very cool calculation). HAPLOVIEW has been used for LD, these authors used it here, what surprised me was they avoided using Arlequin for the LD but used HAPLOVIEW. There's a plethora of recent LD algorithms, never used them so can' say, e.g. quickLD, stuff in R.
What I am definitely saying is the LD works on from raw alleles Thats the most important part of this response. I've attempted to describe the rationale below.

If this is human genomes, or large genome data I'd would recommend plink analysis performed directly on your data, as described here https://www.cog-genomics.org/plink/1.9/ld
Using any sort of matrix will remove the linkage disequilibrium signal, IF its a single matrix for all your data. If its multiple matrices you can perform a correlation statistic between the matrices - but its better just to do a straight linkage disequilibrium calculation, e.g. plink.
The reason why correlation matrices will fail is because LD needs to identify incongruence/disparity and whether that incongruence occurs within random expectations.
Basically you are making a phylogenetic tree - thats what the correlation matrix is about, but the tree assumes clonality. LD is about deviation from random expectations (equilibrium), so either a random distribution, e.g. Poison or a randomisation process is needed. What I'm attempting to say is "a tree" (single matrix) and LD are very different calculations. Whilst, technically you can do LD via matrices (plural), its better to do a LD calculation from the raw alleles.

From the comments, it appears a "tree like structure" (ancestors) is wanted to direct an LD calculation.
What you can do is subdivide the populations and perform LD on each subdivision. A tree (here using a correlation or kinship matrix) is a continuum so it's at what point do you call the the subdivision? Its doable but the information is not presented here.
What you can do is use "Structure" - this has been upgraded alot. This is a HWE, I think using Bayes and it designates how many populations (at LD) are present in your data. You can then perform a specific LD directly on each subdivision.
If you compare the output of "Structure",basically Bayesian HWE you should see a correlation with your kinship matrix. What you are inferring is a partial LD parameter that will define ancestry and whilst its doable it has never been calculated in this manner that I'm aware of: its random or it ain't and thats the criteria of significance. This is definitely how LD works in small genomes (everything below higher eukaryotes).

I think there is a misunderstanding -- I have the raw alleles, and I use them for the calculation, the issue being, since I know for a fact I have families in the data (I just do not know the actual pedigree), I could have a much greater power if this information, or a proxy, i.e. a kinship matrix could be included in the information, because to some degree one could theoretically count how the alleles are passed from (putative) parents to (putative) offspring

@user1256923 answered above