text:
Journal article reference for Illumina's "near monopoly"
Apologies if the question is out of the scope of the bioinformatics SE. Please feel free to close it if it is the case.
Recently, I have been going into a rabbit hole of papers that discuss the progress of sequencing technologies, from the earliest attempts to Sanger, massively parallel sequencing, and long-read sequencing.
To my surprise, I was not able to find any recent reference that direclty highlights Illumina's market dominance. The closest I could find was the 2010 SRA paper ("Almost 80% [...]") that briefly describes the proportion of records per platform or this paper that alludes to Illumina's monopoly. How come it seems so difficult to find recent reports on the use of each sequencing platform?
Am I missing something there, or am I just looking for something that is not that interesting? To be fair, it's probably easier to query the SRA directly to get my answer.

Welcome to the site. Personally (just my opinion) this is much more suited for Biology SE, however someone here may know the answer.

I’m voting to close this question because it is neither really a bioinformatic question, not suited for biology SE. I recommend a google search for "Illumina sequencing market share" which for me yielded some interesting articles, mainly from business-related websites and newspapers. Generally they say that in the (short-read) sequencing field Illumina has a 70-80% share right now which is roughly what I would have intuitively anticipated.

Yes, I'm afraid this isn't on topic here as it isn't about bioinformatics as a science but about the market dominance of a specific company. For what it's worth, this isn't the kind of thing you would expect to find a scientific paper on, if you want a source to cite, look at newspapers, or market analyses. There is no biological point to be made with respect to which company dominates a given market.


text:
What are the output files of RNA-Seq from facility?
This question was also asked on Reddit
I am new in our lab and I am going to do bulk RNA-Seq. What type of files will we get from the company (Genewiz)? Will it be a bunch of Fastq files? or they give something like bam files?

Welcome to the site. Genewiz is a commercial company which you are outsourcing to and I strongly recommend contacting them directly. Essentially you are paying for this service so I think it is fair to request this information directly. Their contact details are here https://www.azenta.com/contact-us ... The implications of the question is that this site provides  technical support for a commercial entity.

Please do not cross-post https://www.biostars.org/p/9571387/ and if you do at least leave a comment with the cross-reference. It's just double effort for people in both communities, especially for these sorts of very basic questions.

OP has deleted their post on biostars.

This cannot be predicted unless you have interacted with them previously. Most commonly a sequencing facility will send raw data as gzipped fastq files, often together with md5 checksums so the user can check file integrity after transfer/download from them.
We worked with a facility before who sent data in unaligned bam format, which then for typical workflows first had to be converted back to fastq.
After all, you're a paying customer so feel free to shoot them an email and ask whatever question comes up. Ask what exactly they will provide so that you can prepare your workflows accordingly. Most likely it will be fastq.gz.

source: Stack Exchange

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

sort: votes


text:
ncbi fastq dump error in loop
I have a very basic objective which I want to give list of id to a ncbi-srafastq tool kit which will prefetch the id which is in .sra file extension and then it will run fastq-dump on the file to split the files into fastq files, single file if there is single end file or paired end file if its paired end or two set of files.
This is my code
#!/bin/bash

input=$1

a=$(basename "$input")
id=${a/.id/}

# Download SRA data
prefetch "$id"

# Extract SRA ID from the downloaded data
#sraID=$(ls -1 | grep -v run.sh)
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')
# Validate the downloaded SRA data
vdb-validate "$sraID"

# Dump fastq data from the SRA file
fastq-dump "$sraID"  --gzip

# Concatenate fastq files for read1, read2, and unpaired reads (if available)
for s in "$sraID"
do
  echo "$s"
  if [ "$s" != "$id" ]
  then
    [ -f "${s}_1.fastq" ] && cat "${s}_1.fastq" >> "${id}_R1.fastq"
    [ -f "${s}_2.fastq" ] && cat "${s}_2.fastq" >> "${id}_R2.fastq"
    [ -f "${s}.fastq" ] && cat "${s}.fastq" >> "${id}.fastq"
  fi
done

# Compress fastq files using pigz
#for i in "${id}"*.fastq
#do
  #echo "$i"
 # pigz "$i" & sleep 1
#done
#wait

# Remove the original uncompressed fastq file if both read1 and read2 files are present
[ -f "${id}_R1.fastq.gz" ] && [ -f "${id}_R2.fastq.gz" ] && [ -f "${id}.fastq.gz" ] && rm "${id}.fastq.gz"

# Test deleting intermediate folder
echo "Space before tmp folder delete"
echo "Removing tmp folder $sraID "
ls *
rm -rf "$sraID"*
ls *
echo 'done'

I run using as small example file as such
./run.sh SRR12332899

The error I get is this
2023-08-05T12:37:32 prefetch.3.0.6: Current preference is set to retrieve SRA Normalized Format files with full base quality scores.
2023-08-05T12:37:33 prefetch.3.0.6: 1) Downloading 'SRR12332899'...
2023-08-05T12:37:33 prefetch.3.0.6: SRA Normalized Format file is being retrieved, if this is different from your preference, it may be due to current file availability.
2023-08-05T12:37:33 prefetch.3.0.6:  Downloading via HTTPS...
2023-08-05T12:40:48 prefetch.3.0.6:  HTTPS download succeed
2023-08-05T12:40:48 prefetch.3.0.6:  'SRR12332899' is valid
2023-08-05T12:40:48 prefetch.3.0.6: 1) 'SRR12332899' was downloaded successfully
2023-08-05T12:40:48 prefetch.3.0.6: 'SRR12332899' has 0 unresolved dependencies
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Table 'SEQUENCE' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'ALTREAD': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'QUALITY': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'READ': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'SPOT_GROUP': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'sra/SRR12332899.sra' contains only unaligned reads
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' is consistent
2023-08-05T12:40:49 fastq-dump.3.0.6 err: item not found while constructing within virtual database module - the path 'sra' cannot be opened as database or table
fastq-dump quit with error code 3
sra
Space before tmp folder delete
Removing tmp folder sra
run.sh

sra:
SRR12332899.sra
run.sh
done

If I do just
prefetch SRR12332899 and then fastq-dump SRR12332899
I get a SRR12332899.fastq  SRR12332899.sra
This is for one file since I want to download list of files using the script which is not working.
Any suggestion would be really helpful

Note that you are catting an already gzipped files since you used --gzip in fastq-dump, so the suffix of the output must be .gz
It would be muchfaster to simply search for the dataset at sra-explorer.info and get fastq download links right away.

The main problem is that your script is inadvertently treating a set of file and directory names as a single SRA accession number, all in one variable.  You should update it to work with each accession separately and make sure they really are accessions.

On this line you're taking all the names of things in your working directory (minus the script), each on a separate line, and trimming .sra off the end of each line:
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')

In your case that list was apparently just a directory named "sra".  (Even though I don't have that directory on my end, you can see it in your output, from the first of two ls * at the end, where you have sra/SRR12332899.sra.)  So later on it tries to fastq-dump sra, which gives you that "cannot be opened as database or table" error.  (vdb-validate can take a directory as input, so it doesn't complain.)  If I run your script from scratch in an empty directory with SRR12332899, it actually seems to run OK, aside from removing everything it creates at the end.  But if you run it with a second accession you see the problem, since it then tries things like vdb-validate "several accessions and filenames separated by newlines".
I think a proper fix would involve a bunch of changes, but here are some tips, at least:

use subdirectories to keep things organized so you don't have to juggle multiple types of file in one place
set -x at the start of a script will show you each command as it happens, which can be very helpful for troubleshooting
be careful with quotes and whitespace in bash.  If you want to loop over each word in a variable you can do for thing in $var; do ...; done (unquoted) but if you run into filenames with whitespace, it breaks them up.  Another pattern is command | while read thing; do ...; done which will go line-by-line.  (In that case if you were working with a variable storing the command output, you might do something like echo "$output" | while..., quoted that time, since you'd want to preserve the newlines.)
ls * will expand each file/directory name as separate arguments (so bash really ends up calling ls "sra" ...whatever else...), and when ls sees directory names is lists their contents too.  If you just want to show all the files/directories you can leave off the asterisk altogether, or if you wanted to give directory names but not have it list their contents, you can use ls -d.
you don't necessarily have to call ls to store the names as a block of text.  You can loop over each filename, including those with spaces, directly in bash (e.g. for thing in *; do echo "path: '$thing'"; done)
be careful with the rm command, especially combined with wildcards, and especially with -rf.  If sraID happened to be empty somehow, think of what rm -rf "$sraID"* would actually mean.  (And yes, scripts can remove themselves!)

will work on your suggestion and show you what i have done

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


text:
Interpretation of clusterProfiler results
I am using clusterProfiler for differential gene expression using this code:
https://yulab-smu.top/biomedical-knowledge-mining-book/enrichplot.html
I have a question about the third line of code:
 de <- names(geneList)[abs(geneList) > 2]

Here all gene names are taken with absolute fold change >2. As an example, if in a hypothetical extreme case you have 100 downregulated genes (log2 fold change <-2) and zero upregulated genes (log2 fold change >2), then you may incorrectly think that identified pathways are activated, but they are actually repressed. How is this accounted for? Thank you.

Hypothetical scenarios don't make the best questions, but this question's ok. Do keep in mind this is differential gene expression rather than explicitly being repression (don't always intersect).

I agree with you :)

The scenario you described would be reported by the code because:
... [abs(geneList) > 2]

The abs will change any negative value to positive and retain all positive values as positive. Therefore values <-2 will be reported as well as values >2. abs = absolute and to take the absolute value of a number is exactly the operation just described.

From comments. Just remove abs
negative <- names(geneList)[geneList < -2]
positive <- names(geneList)[geneList > 2]

Thank you for your feedback M__. Yes, this operation selects all genes with >2 and all genes with <-2 fold change in expression. However, none of the values are themselves retained. It is only the names of the genes that are taken forward to the next step. Therefore, the direction (- or +) is lost after this step. In other words, if you have gene A +8 fold change, gene B +3 and gene C +2.4, then the end result is the same if gene A had been -6.2 fold change, gene B -2.4 and gene C +2.1. I suppose I conclude that we are looking at representation, rather than at up/downregulated pathways?

@SylviaRodriguez addressed in the comments

Thank you. I like that solution. Two separate analyses.

Thanks @SylviaRodriguez I'm glad its resolved.


text:
spliceAI not giving output value while running using vep (Variant Ensemble Predictor)
While generating the predictions using spliceAI from vep, I am not getting any values.
I used this command to generate the output,
vep \
  -i input.vcf \
  -o output.vcf \
  --cache \
  --force_overwrite \
  --tab \
  --no_check_variants_order \
  --force \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz

Below, are the examples of entries in the output file.
ENSEMBL VARIANT EFFECT PREDICTOR v109.3
Output produced at 2023-08-07 09:05:43
Connected to homo_sapiens_core_109_38 on ensembldb.ensembl.org
Using cache in /.vep/homo_sapiens/109_GRCh38
Using API version 109, DB version 109
ensembl-io version 109.4946a86
ensembl version 109.10baaec
ensembl-funcgen version 109.cba2db8
ensembl-variation version 109.18a12b6
...
Uploaded_variation Location Allele Gene Feature Feature_type Consequence cDNA_position CDS_position Protein_position Amino_acids Codons Existing_variation IMPACT DISTANCE STRAND FLAGS SpliceAI_pred
100017682 chr1:100017682 A ENSG00000117620 ENST00000370153 Transcript missense_variant,splice_region_variant 981 880 294 A/T Gca/Aca - MODERATE - -
100017682 chr1:100017682 A ENSG00000117620 ENST00000370155 Transcript splice_region_variant,3_prime_UTR_variant,NMD_transcript_variant 882 - - - - LOW - 1 - -

Is there any wrong in the input file? Could anyone help here?

Thanks (again) for editing this @RamRs. We need to make documentation (meta-post) for new member's code format .... in fact this is particularly well edited and like to use this as the example for a meta post.

I use a template at biostars: https://www.biostars.org/p/9571077/#9571078


text:
Using CSV files continaing scRNA-seq count data from GEO
I'm completely new to scRNA-seq analysis (and to most of Bioinformatics as well) so I apologize if this was asked a million times before.
I am trying to get the single cell expression data of lymph node and bone marrow from this GEO entry
In Samples, there are multiple organs and accessing any of them gets me a CSV file compressed in gz format. The CSV file appears to contain the barcodes and the counts, but I do not know where to go from there. I need to annotate each cell type and I am not sure how. Seurat tutorials indicate that I have to have a few other file types which do not seem to be available in this repository. Does this mean I cannot use it, or is there a way to get this data to count/cell type just using what's available? I imagine that I don't need to perform any alignment since the rows already contain gene names.
Any and all guidance would be very much appreciated. Thank you so much in advance.

Welcome to the site @AhmadAl-Shaibi. My advice is to read the Seurat tutorial, attempt the calculation on your test data set and then report back with specific errors.

Please do consider revising the question. Just to mention there are close requests for this question.

Please understand that your question is open-ended and essentially you're asking for a complete guidance through your analysis which is not feasible. I recommend to go through either the Seurat tutorials or the Bioconductor single-cell workflow to learn the basics, such as typical workflows, formats of required input files etc.
With this knowledge you can then start finding out which sort of data are provided in GEO and start getting your analysis going. You're welcome to ask specific questions here, but please be precise and not open-ended.


text:
How to find the file names referred to in the R script produced by GEO2R in GEO Accession
I am new to both coding and GEO so this may be a silly question.
I need the raw_counts data set for the GEO data set I am using but data in the downloads box is not the raw_counts. How would I be able to find the raw counts?
As I am new to both I decided to use GEO2R to run an analysis and then to see if I could create the same results using code in R as I new that I could look at the R script for GEO2R if needed.
I was struggling with which of the files I needed to download so I looked at the R script used in GEO2R but the files that are being used are completely different to the ones that I have downloaded.
I am using the data set for Series GSE178646 "CRISPR gene editing and inducible pluripotent stem cell neuronal disease modelling for rare disease diagnosis: EMHM1 genetic variant analysis in Kleefstra Syndrome" and am only interested in CF7-12 but the files that I have downloaded are called "GSM5395222_CF7_kallisto_abundance.tsv.gz" where as in the R script in GEO2R the code has used the file "GSE178646_raw_counts_GRCh38.p13_NCBI.tsv.gz" and I am really not sure where I could access this file.

Apologies if there appear to be more than one question the only question is how to find the raw counts data in for this dataset. How would you suggest editing the question so it is clear thats the only question

Thanks for the edit @Sabrina, your edit is clearer. The post is fine - it comprised a singular question when I went through it - and was my mistake, apologies.

All the data associated with this project is here:
https://www.ncbi.nlm.nih.gov/gds/?term=GSE178646
The actual download page and project description is here
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE178646
The dataset you seek GSE178646_raw_counts is here:
https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE178646&format=file
The NPC EHMT1_SNV CF7 cells for GSM5395222_CF7_kallisto_abundance.tsv.gz is on the same page here.
If you now click on the GEO2 analysis on this page it will take your here where you select your samples and begin the GEO2 analysis.

From the comments,
How to do this?
Its easy simply go to https://www.ncbi.nlm.nih.gov
In the check box search for the option GEO datasets, in this example I queried GSM5395222 in the box. This goes to the search results here, the top hit called "1" is document you want here. The very top hit (above the 1) is actual data in question, minus the project information page.
All the information is within this document and simply use a "search text in webpage" using GSM5395222

Thank you for this, how were you able to find the raw counts data or do I just put raw counts and it will find the data?

@Sabrina I've modified the post above to include this information below the --------

Thank you very much for your help


text:
Multiple SRR ID downloaded for single GSM IDs input issue
So here I'm trying to donwload multiple fastq files based on the GSM ID input the script works fine when then input is only SRR id but it runs into error when it is given GSM ID
My small GSM ID input set
cat srr_id.txt
GSM5369701
GSM4663476

I only see two GSM IDs in your output, is it the extra SRR* files you don't want? I don't really have much experience with the SRA, but I can see that, for example, SRR12180300 is under GSM4663476: https://www.ncbi.nlm.nih.gov/sra/?term=SRR12180300. Are you sure there is an issue here?

yes I mean if two GSM ID then 2 SRR file I should be getting but which is SRR12180300 and SRR14836042 for the other GSM ID but Im getting other also why?

@terdon these are technical reads which I don't want in other words i need final output where I get only one output instead of different technical reads i thought fasterq-dump would do this by default but it it not doing it

Basically you're asking why a given GEO accession (GSM...) has many Runs (SRR...). That is because the submitters uploaded multiple fastq files for the same sample, which is the case if a library has been sequenced over multiple lanes or flowcells, often to increase depth. You can use cat to combine these "technical replicates" into a single fastq file for both R1 and R2 file in case of paired-end data, or one single fastq file in case of single-end data.
For example, if you have file1.fq.gz, file2.fq.gz and file3.fq.gz just do:
cat file1.fq.gz file2.fq.gz file3.fq.gz > combined.fq.gz

"You can use cat to combine these "technical replicates" into a single fastq file for both R1 and R2 file in case of paired-end data". it is okay to go head with this approach meaning it wont be an issue with the downstream alignment .?

cat is fine, see edit.


text:
Bioprocess data visualisation and analysis
I'm looking for some kind of software, ideally FOSS, that will allow us to plot, compare, annotate and add additional points to time series data from bioprocessing runs (bioreactors).
Many of the bioreactors come with some of their own quite limited software to do this, but it cannot be installed on other computers.
I have looked around a fair amount but can only expensive-looking biopharma software. Is there any academic software to achieve this?
I was considering that writing it myself wouldn't be that difficult with python, but I don't want to waste effort.
(maybe this SE is not the best for bioprocessing, not sure where else to ask.)

Welcome @Alext. Generally, the question isn't sufficiently specific about what is required, few of us will understand bioprocesser measurements, although alot will understand timeseries. In the public domain one approach is adapting IoT https://www.mdpi.com/2076-3417/11/24/11932, where IoT is connecting a device to a computer and accessing the output. GCP has strong IoT capabilities.

As a general answer the technique here is called "Internet of Things" (IoT), where a machine, or an "IoT" device reports to either the cloud or another computer. In bioinformatics this technology has not embraced at present outside the industrial sector.
Google Cloud Platform (GCP) does have an IoT platform where the devices can be registered and data is streamed to the cloud. However, doing this operation across a LAN (local area network) is automatic and can be directed via cron notably in automating and collecting the data and analysing it periodically. For a generic IoT stack consider https://www.mdpi.com/2076-3417/11/24/11932

Thanks for your response :), however I am not asking about receiving the data in real time in this instance, I am aware of some industrial standards like OPCUA that can be used for this, but rather some software to use after the run, to view and analyse the recorded data.


text:
phytools add.species.to.genus could not match your species to a genus
I have a species tree to which I am trying to add species, which includes a tip "Geissospermum laeve". When I run:
with_polytomes = phytools::add.species.to.genus(tree, "Geissospermum polytome1", where='root'), the original tree is returned with the error message: 'could not match your species to a genus  check spelling, including case'.
Reproducible example:
test_tree = ape::rtree(5)
test_tree$tip.label = c("t3 sp1", "t4 sp1", "t5 sp1", "t2 sp1", "t1 sp1")
with_polytomes = phytools::add.species.to.genus(test_tree, "t1 polytome1", where='root')

I don't understand why this is happening, if anyone could shed some light on what's causing this error that would be great

Two things to try,

setting genus=NULL, The logic being it doesn't know which genus to add it to. Just as a precaution setting all blank space within a taxa id in the tree to _ ...

add.species.to.genus(test_tree, "t1_polytome1", genus=NULL, where=c("root"))

Maybe replace root with random to get it working.

alternative using tree instead

genus.to.species.tree(test_tree, "t1_polytome1")

I tried lots of different combinations of parameters but in the end the only thing that resolved the issue was using underscores instead of spaces in the tree labels i.e. test_tree$tip.label = c("t3_sp1", "t4_sp1", "t5_sp1", "t2_sp1", "t1_sp1"). Curiously though when the tip labels are set like this, the function works  with species = "t1_polytome1" and species = "t1 polytome1".
WRT (2), phytools::genus.to.species.tree(test_tree, "t1_polytome1") returned a tree without t1_polytome1 added and no error given, maybe because it expects  a genus-level tree.

Thanks @A.Bollans thats good debugging. In tree files removing all blank space is standard practice, in phylip format this a bug.


text:
Help with finding an unsupervised learning algorithm for RNA sequence/structure pattern identification
I have a bunch of RNA sequences (and their optimal secondary structures) and their corresponding energy values (measured by mean free energy) and I'm trying to find a way to identify features (patterns in their structures or sequences) common between samples of similar energy values. the overall goal is to try and identify sequence or structural elements in the RNA samples that tend to correlate with low mean free energy values. would be grateful if anyone could recommend me algorithms to look into using - I'm assuming this would be an unsupervised learning project and I only have experience with supervised stuff, but I'm looking into PCA right now and not sure if that'd be useful. should I be looking into something else?
Intuitively I'm imagining it as like a clustering + feature extraction problem where I have a bunch of dots, each representing an RNA sample, and then the axis represents energy so the dots could be clustered in energy value similarity and then within each energy value cluster patterns/relationships could be found between the sequences and structures of the samples. but not sure if an ML algo exists to handle this and PCA seems like not what I'm looking for because the axes would be principal components and not energy...

It depends what is meant by "features". Features can be a particular structure component in the RNA structure. In ML/DL these "features" are called instances (or observations). "Features" is also a technical word used in machine/deep learning, which can be completely different, it means the data or statistical variables associated with an instance.
Feature weight extraction (ML/DL definition) using unsupervised learning is complicated, but doable. This is where the features corresponding to each entry are weighted, i.e. quantified. Under traditional unsupervised learning this is not achievable, unless you want to include/exclude "features" and assess cluster integrity. I personally considered it a lot of work for no real interpretable results.
If the features are components of the RNA structure (what would be technically called 'instances' in ML/DL) then the general traditional unsupervised learning following PCA analysis is tSNE. However, tSNE is often used to get better resolution the PCA feature reduction output, i.e. tightens the clusters. PCA->tSNE is a cool technique and works well. PCA and tSNE can give different clusters by the way. As a rule of thumb I consider PCA > tSNE.
Multi-dimensional scaling etc ... needs a defined analytics aim.

Information on tSNE and PCA this online resource looks fine
https://www.datacamp.com/tutorial/introduction-t-sne
The above is if you want to do tSNE and compare the results to PCA this link is good.
In terms of using tSNE to sharpen PCA clusters I always used O'Reilly's "Introduction to Machine Learning with Python". There is a very clear example of Python code of the technique there, towards the end of the book and why it works.

Thanks! Do you have any resources I could look into that would help me understand how to frame this as a PCA problem as someone new to it? If not no worries

Makes sense re: edit, thank you!

Dear @joshua213p I have answered below the --------

Thanks @joshua213p

hey! reading through the t-sne blogpost and I'm a little unclear on how i'd use either of these algorithms here for my goals. i get it'd cluster my samples based on similarity in terms of first and second principal components/t-SNEs, but what could i then use identify structural or sequence patterns within each cluster? please let me know if I'm missing something - thanks!

Hi @joshua213p can you submit this as a separate question please?


text:
detecting DNA of possible infection on cell line
we've done Whole Genome Sequencing with Illumina on DNA extracted from cell lines, where some are suspected to have been infected by something (viral, bacterial, not sure).
I've tried running metaphlan4 (version from github) and got zero hits in all of them.
What other tools could I be trying? Thx.

I've had some success looking for contamination in Illumina and Nanopore reads using kraken2.
Kraken2 has a phylogenetic-based approach, which in the first instance tries to find the most specific taxa that the target read likely comes from (going further up in the tree if there's ambiguity), and also produces a summary table indicating the approximate proportion of reads that come from particular taxa. It can be paired with bracken for additional reliability.

Yes I 100% support this approach. Kraken is an excellent algorithm and will identify a whole diversity of contaminants. There are contamination approaches IF the contaminant is predicted. When its anything Kraken is an extremely good approach. Setting up the Kraken database is a pain, but this case it should be easy because you doesn't have to be very accurate , you don't need the species of the contaminant. Absolutely commend this approach. I don't think bracken is needed in this particular case - outside this case it is needed to correct the count distribution bias.

Maybe less advanced but worth mentioning, I've had good results with fastq_screen. It's simply trying to align a subset of your sequences against various genomes (comes with a few genomes, you can add more).

I didn't manage to try the ones reported by @gringer and @alexlok here, but I did manage to perform an analysis using the OneCodex webapp: https://app.onecodex.com/ and based on the experimental information for the samples being analysed, the results were concordant with the expectation. You can do quite a lot on the "free tier", so at the moment I didn't upgrade to the paid version.


text:
How to get coefficient tables from multiple regression model result!
I've done multiple regression modeling in R using the below script. Now I'm wondering how can I make a table for the coefficients of each model of all models from lm, glm, and rf separately.

all_variables <- names(PaC09c_w20_dat3)
response_variables <- all_variables[c(1:3)]
predictors <- all_variables[-c(1:3)]

lm_model <- lapply(
  response_variables,
  function(x) lm(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3)
) |>
  setNames(response_variables)

glm_model <- lapply(
  response_variables,
  function(x) glm(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3)
) |>
  setNames(response_variables)

rf_ranger <- lapply(
  response_variables,
  function(x) ranger(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3, importance="impurity")
) |>
  setNames(response_variables)

save.image(file="PaC09c_w20_Reg.RData")

Input data looks like (Original Table is very large),
  head(PaC09c_w20_dat3)
        EE87865ln1   EE87866ln1   EE87895ln1  blood_vessel_w20 adrenal_gland_w20 bone_element_w20 bronchus_w20
               <dbl>       <dbl>      <dbl>     <dbl>             <dbl>            <dbl>        <dbl>
        1    0.00391     0.00326    0.00332       0                 0                1            0
        2    0.00139     0.00116    0.00132       0                 0                0            0
        3    0.00360     0.00270    0.00469       1                 1                0            1
        4    0.00323     0.00348    0.00339       0                 0                1            0
        5    0.00323     0.00330    0.00382       0                 1                0            0
        6    0.00278     0.00208    0.00214       0                 0                1            0

The expected output would look like this,
                           EE87865_lmcof   EE87866_lmcof
blood_vessel_w20           6.332633e-06    5.916615e-06
adrenal_gland_w20         -2.159470e-06   -2.399176e-06
bone_element_w20           8.578804e-06    8.106153e-06
bronchus_w20              -2.956291e-05   -2.152231e-05

Your code is neither reproducible (one would need the bed file) nor small enough to mentally compute what could be in each variable. How do you expect people on the forum to help you get to a random goal you have?

Well, I found a solution,
 coefs_lm<-as.data.frame(lapply(lm_model,function(x)coef(x)[2:5]))

 coefs_glm<-as.data.frame(lapply(glm_model,function(x)coef(x)[2:5]))

var_imp <- function(object, ...) {
  var_imp_r <- object$variable.importance
  var_imp_r
}

coefs_rf<-as.data.frame(lapply(rf_ranger,function(x)var_imp(x)[1:4]))

2:5 for lm and glm is the number of samples EE.., where 1 is intercept.

What is the 2:28 part? Is it specific to your dataset(s)?

Sorry, it should be 2:5 for this toy dataset in this example.

To get the coefficients alone
lm_model$coefficients
glm_model$coefficients

Everything including coefficients its,
summary(lm_model)
summary(glm_model)

I assume this is the same for summary(rf_ranger), albeit not sure whether coefficients would work here

lm_model$coefficients and glm_model$coefficients produces NULL. As you can see every model object (lm_model, glm_model, and rf_ranger) consists of multiple regression models of samples (e.g., EE...).

Given how OP has made these things lists, we'll probably need an lapply there.


text:
Applying PCA or tSNE on non-numerical data
Reading through this t-sne blogpost and I'm a little unclear on how i'd use either of these algorithms here for my goals. i get it'd cluster my samples based on similarity in terms of first and second principal components/t-SNEs, but what could i then use identify structural or sequence patterns within each cluster? please let me know if I'm missing something - thanks!
I'm also curious as to how I would apply PCA or t-SNE on my data considering its non-numerical nature. I was thinking of applying one-hot encoding, but is there a generally more suggested approach? I've attached a screenshot of a few of my dataset rows for context (for the sake of having the table fit, I chose rows with simpler structures/zero MFE values, but most rows have complex structures and non-zero MFE values). Thanks.

This is based on my previous post:

I have a bunch of RNA sequences (and their optimal secondary structures) and their corresponding energy values (measured by mean free energy) and I'm trying to find a way to identify features (patterns in their structures or sequences) common between samples of similar energy values. the overall goal is to try and identify sequence or structural elements in the RNA samples that tend to correlate with low mean free energy values. would be grateful if anyone could recommend me algorithms to look into using - I'm assuming this would be an unsupervised learning project and I only have experience with supervised stuff, but I'm looking into PCA right now and not sure if that'd be useful. should I be looking into something else?

Intuitively I'm imagining it as like a clustering + feature extraction problem where I have a bunch of dots, each representing an RNA sample, and then the axis represents energy so the dots could be clustered in energy value similarity and then within each energy value cluster patterns/relationships could be found between the sequences and structures of the samples. but not sure if an ML algo exists to handle this and PCA seems like not what I'm looking for because the axes would be principal components and not energy...

I'd recommend PCA, tSNE the PCA (first)+ tSNE (second).
The sequence data - one hot key encoding is a good idea.
import pandas as pd
datadict = {'UTR':[.....],'Sequence':['AGAA','ACAG' ...], 'Optimal Secondary Structure':[....]}
df = pd.DataFrame(datadict)
df_onehot = pd.get_dummies(df['Sequence'])
df = pd.concat([df, df_onehot], axis=1)     # add the one hot encoded column to master dataframe
df.drop(...) # remove the original 'Sequence' column

The other issue is vectorising the "dots". This is really a separate question, here's what I'd do ...
For optimal secondary structure, maybe simply count the dots?
count = df_series.str.count(".")

This only works on a series so you'll need to do something like df_list = df['Dots column'].to_list() and then df_series = pd.Series(df_list). Finally the count will need to be concat back onto the df and the original Dots deleted.
The thing I'm not clear about is your 'instances' that could be 'UTR_info' thats the stuff you are calling you data or "features".

Got it, but how does pca and tsne get me closer to actually identifying structural/sequence patterns present in the created clusters?

Thanks. You can now upvote BTW. This is about the fundamental design of your experiment. If post another question I'll sketch an outline about what unsupervised learning achieves. I don't exactly know your project aims (I can guess) but I know what unsupervised learning delivers.


text:
What algorithm should I use to find what RNA sequence/structural elements lead to different labelled MFE values?
To put it super simply: I have RNA sequences, their optimal secondary structures, and the MFE values of those structures. I'm trying to computationally identity what structural elements and/or sequence motifs tend to be present in RNAs with lower MFEs.
Was thinking about using unsupervised learning for this, but is this the wrong approach?
Past version of this question

I would maintain traditional unsupervised learning is a very good method here (PCA and/or tSNE). The number 1 question is which column is the observation (label)? That depends on the aim of the study. The absence of a "target" column is the thing that separate unsupervised learning from supervised learning.
The "observation" (label), i.e. the column that provides the name, BUT doesn't take part in the calculation, will either be a gene or an isolate/species/taxonomic-group and the aim is to see which genes or isolates cluster together. The implication is that the MFE (MLE?) is the observation (label) ... can be I guess. Again its about the aim of the study.
The trick with unsupervised learning is to get as many "features" (columns of data) for each "observation" (probably gene name or taxonomic group). MLE (MFE?) energy (entropy?) has a value, but there are other values (which I forget but other entropy values), those can be included. Dumbie variables are used for e.g. sequence data, but if there are counts will work and experimental data measurements can be included. There are a few rules:

When you have a cateogorical variable as a column in a dataframe it must be designated as such, pandas will allow this. Similarly if you have continuous variable, MLE that must be designated as a float and counts will be designated as integers. You don't want a float to appear as a categorical variable, that will be disastrous. You need to set the data type of each column (see below).
In PCA a reasonable number of components must comprise a large amount (sum) of the variance, e.g. >70%.
You can only see realistically two components at any one time on a graph, just plot them out, e.g. PC1 versus PC2, then PC1 vs PC3. Look at which give nice clusters. Don't just look at PC1 and PC2 and believe thats the only answer (which is what everyone does).

Point 1 is important if the PCA or tSNE may not permit mixed data types (I can't remember), but if so it will complain.
Point 2 If the sum of the first three principle components exceed 70% of the variance - thats very, very good. If its alot of then you might  consider performing tSNE on the PCs to tighten the result.
Thats about it, what you are looking for is clusters that give some indication about the way the data is behaving. This can be the end point of the analysis or it can be the hypothesis for further testing.
Very clear distinct cluster which make biological sense is a good result with unsupervised learning. Obviously the issue is that the clusters are not significant, but if they look good and make sense that's the key to the analysis.


text:
Importing local files to a GATK docker image
I have pulled GATK docker on my computer but when running a GATK command, can not locate my input file
My OS version
cat /etc/os-release
NAME="Ubuntu"
VERSION="18.04.6 LTS (Bionic Beaver)"
ID=ubuntu
ID_LIKE=debian
PRETTY_NAME="Ubuntu 18.04.6 LTS"
VERSION_ID="18.04"
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
VERSION_CODENAME=bionic
UBUNTU_CODENAME=bionic

Files in my working directory
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ ls
GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDup.bam
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDuplicates.txt
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDup.bam
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDuplicates.txt
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDup.bam
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDuplicates.txt
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$

I then pulled GATK docker
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v ~/vcf:/gatk/vcf -it broadinstitute/gatk:latest
(gatk) root@8b2af4d7aadf:/gatk#

(gatk) root@8b2af4d7aadf:/gatk# ls
GATKConfig.EXAMPLE.properties  gatk-package-4.4.0.0-local.jar  gatkPythonPackageArchive.zip  scripts
README.md              gatk-package-4.4.0.0-spark.jar  gatkcondaenv.yml          vcf
gatk                   gatk-spark.jar              gatkenv.rc
gatk-completion.sh         gatk.jar                run_unit_tests.sh
(gatk) root@8b2af4d7aadf:/gatk#

As you see, folder seems to be empty because I can not connect my local working directory to container
(gatk) root@6d3d3ff69e21:/gatk# cd vcf
(gatk) root@6d3d3ff69e21:/gatk/vcf# ls
(gatk) root@6d3d3ff69e21:/gatk/vcf#

I have tried all of these but did not work
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v /data/Continuum/WES/vc:/gatk/ -it broadinstitute/gatk:latest

fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v ~/data/Continuum/WES/vc:/gatk/ -it broadinstitute/gatk:latest

What does -v ~/vcf:/gatk/vcf do? Why are you using ~/vcf there? ~ is your $HOME directory, is /data/Continuum/WES/ your home directory? Does it work if you try docker run -v /data/Continuum/WES/:/gatk/vcf  -it broadinstitute/gatk:late?

Yes you solved that, thanks a lot

@terdon answered that
docker run -v /data/Continuum/WES/:/gatk/vcf -it broadinstitute/gatk:latest

Glad it's solved. Unless @terdon posts this as an answer, could you kindly mark this as "resolved" please (this helps with site stats). The "resolved" flag can always be shifted between posts, i.e. if terdon answers later.

I do not know how to do that sorry

Just "accept" your answer (see https://bioinformatics.stackexchange.com/help/someone-answers) by clicking on the checkmark on the left.


text:
Interpreting GFA graph visualized in Bandage
I assembled Nanopore sequenced reads with Flye and visualized the GFA graph in Bandage but I don't really know how to interpret the result. For context, this is a yeast (DNA) genome. My goal is to figure out areas of repeats, ploidy level, and overall structure of the strain. So far I understand that the bubbles are repetitive areas but I mainly want to know how to interpret a graph like this.

I'm new to interpreting these graphs and I couldn't find any sort of guide online so if someone has any resources they'd be willing to share or if they could help interpret this graph I would really appreciate it!

Please [edit] your question and add some detail. To be honest, I don't have a clue what this is, I've never seen this sort of thing before, so your question might be clearer to someone who knows more, but adding some detail can't hurt. For example, what organism did you sequence? What kind of sequencing (RNA? DNA?). Also, why did you plot the GFA? What question are you asking of your data?

I've used Bandage to visualise transcriptome assemblies, but I'm not sure how to answer your question/s about repeats, ploidy, and structure sorry. If you don't get an answer here, I suggest reaching out to the author - Ryan Wick (email: ryan.wick@unimelb.edu.au / mastodon: @rrwick@genomic.social / github: https://github.com/rrwick/). Good luck!

I would suggest making a dot plot of your assembly against a reference (related) yeast genome assembly.
I'd suggest also running QUAST against the yeast reference to get a finer grained idea of the comparison.
Bandage is more of an assembly QC tool for understanding whether there are issues with the assembly. It will not tell you anything about ploidy or that kind of thing.
Off the top of my head, I find it worrisome that 90% of your sequence is in a single scaffold/contig. Presumably there's more than a single chromosome. It does look like the repeats are messing you up.
I'd suggest using hifiasm instead and seeing if it makes different decisions.


text:
Parsing gene names to detect organism from h5ad input
I have one function that is for reading .h5ad files
ReadH5AD <- function(file) {
  # Load libraries
  library(anndata)
  library(reticulate)

  # Check if the file has the .h5ad extension
  if (tolower(tools::file_ext(file)) != "h5ad") {
    stop("Only h5ad file paths are accepted!")
  }

  # Import and read the h5ad file
  ad <- import("anndata", convert = FALSE)
  h5ad <- ad$read_h5ad(file)

  # Extract counts matrix as a dense matrix
  counts <- as.matrix(py_to_r(h5ad$X))

  # Extract cell names
  cell_names <- py_to_r(h5ad$obs$index)

  # Extract gene names
  gene_names <- py_to_r(h5ad$var$index)

  # Transpose the matrix
  counts <- t(counts)

  # Set cell names and gene names as row and column names
  colnames(counts) <- cell_names
  rownames(counts) <- gene_names

  return(counts)
}

Now there is another function which is used for file which are .mtx type along with barcode.tsv and gene.tsv.
In order to detect organism type this function is used
detectOrganism <- function(outdir,data){
  gene.loc <- file.path(outdir, list.files(path=outdir, pattern='*genes\\.tsv', recursive=FALSE)[1])
  features.loc <- file.path(outdir, list.files(path=outdir, pattern='*features\\.tsv\\.gz', recursive=FALSE)[1])
  tryCatch({
    # Flag to indicate if data is from CellRanger >= 3.0
    checkGeneFile <- file.exists(gene.loc)
    if (!checkGeneFile) {
      genes <- read.table(file = features.loc, header = FALSE, sep = '\t',nrow=1, row.names = NULL)
    }
    else{
      genes <- read.table(file = gene.loc, header = FALSE, sep = '\t',nrow=1, row.names=NULL)
    }
    if(grepl("^ENSMUSG",genes[1,1])){data$organism="Mouse"}
    else if(grepl("^ENSG",genes[1,1])){data$organism="Human"}
    else{data$organism="NA";stop()}
  },
  interrupt = function(x){errorOccured<<-TRUE; speciesErrorNotification()},
  error = function(x){errorOccured<<-TRUE; speciesErrorNotification()}
  )
}

This works when the input type is matrix.mtx,gene.tsv along with barcode.tsv. As this can detect using grepl from the gene.tsv file.
I would like to implement the same for .h5ad input.
  gene_names <- py_to_r(h5ad$var$index)

This part can extract the gene_names from a given .h5ad input file but how do I pass this to the detectOrganism function.
I have here basically two file types .h5ad and .mtx.
So my simple goal is the add the same functionality when the input is .h5ad instead of .mtx file type.
Any suggestion or help would be really appreciated

Please [edit] your question and add some example input files that we can use to reproduce the problem. Also show us the expected output from those examples so we can know if we get it right.

do i need to add the h5ad file?


text:
Why do we delete scaffolds shorter than 500 - 1000 bp from the assembled genome?
After assembly of genome, some protocols sometimes call for removal of scaffolds shorter than 500 or 1000 (some papers have one number while the other has the other.) Is this simply to remove the artefact and to make later analysis smoother, without the extroneous data? Also how do we know how where is our cut off point for size?

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

I recognise my answer is late but I hope it helps. The following has no emperical proof, just anecdotal from working with bacterial genomes:
I'm not sure what kind of genome you're working with, but from my work with bacteria, we often delete scaffolds shorter than 1000 - 3000 bp because these are considered "noise". Essentially the scaffolds are too small to gather any particularly useful information from as they may only carry a single gene (or a handful of little genes). If you've added in long read data with your short, half the time those small scaffolds are just artifacts from sequencing errors.

Oh and for the cutoffs: I work with metagenomes, so I tend to work my cutoff around my assembly N50. Say the N50 is 2100, I will round that to 2000 and remove any scaffolds smaller than that. It seems to help me save the majority of the data whilst minimizing noise. Each case will be different though :)

Agree with what @rainman has said. depending on the organism, you might want to go lower, eg if the genome is is tiny you might drop the cutoff. Also it's good to note there are scaffolding tools that help bridge the gap between contigs.


text:
WGCNA error: Length of colors vector not compatible with number of objects in 'order'
plotDendroAndColors(bwnet$dendrograms[[1]], cbind(bwnet$unmergedColors, bwnet$colors),
                    c("unmerged", "merged"),
                    dendroLabels = FALSE,
                    addGuide = TRUE,
                    hang= 0.03,
                    guideHang = 0.05)

I am getting the following error. colrs modules are not coming in  dendrogram graph. I am using 4gb ram laptop. is this problem because of laptop or my mistake in understanding code.
Error in .plotOrderedColorSubplot(order = order, colors = colors,
rowLabels = rowLabels, : Length of colors vector not compatible with
number of objects in 'order'. madam getting this error.

My total genes are 15456. Under blockwisemodule I chose 5000

Welcome to the site. Firstly, what data set is this? Weighted correlation network analysis (WGCNA) whilst published in a bioinformatics journal could be applied to any data set, i.e. not necessarily a biological data set. The context of this question is quite important. WGCNA (2008) is not widely used and there could be alternative approaches.

Please follow the code in WGCNA tutorial section 2c. Since the analysis ran in blockwise manner, you need to plot the dendrograms for each block separately and subset the module colors to each block. The example code is in the tutorial.


text:
Invalid tag name: "1KG" and Invalid character '/' in 'SA' FORMAT field at chr1:6197766
I have a VCF file that is of the following format:
##fileformat=VCFv4.1
##fileDate=20210329-231145
##AnalysisPipelineName=NGS5
##AnalysisPipelineVersion=V5.2.1
##AnalysisConfigurationName=NGS5_Exome
##AnalysisConfigurationVersion=5.2.1
##NGSPanelName=SureSelect_ExomeV7_plus_720G
##NGSPanelVersion=V1.0
##BatchID=WES_032721_A
##Lane=Lane1
##reference=hg38
##PathogenicSilentFile=Pathogenic_Silent.20200730.vcf
##FlowcellName=210327_A01233_0050_AHYFJHDSXY
##VariantDefDateTime=20210323-070701
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=VF,Number=1,Type=Float,Description="Variant Frequency">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allele Depth">
##FORMAT=<ID=SB,Number=1,Type=String,Description="Strand Bias with read numbers (Fw/Rv)">
##FORMAT=<ID=SA,Number=4,Type=Integer,Description="Number of 1) forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles, used in variant calling">
##FORMAT=<ID=SP,Number=1,Type=Float,Description="Phred of Chi-square p-value for strand bias">
##FORMAT=<ID=SO,Number=1,Type=String,Description="abs(log10( (RefFw*AltRv)/(RefRv*AltFw) )) for strand bias">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=TI,Number=.,Type=String,Description="Transcript ID">
##INFO=<ID=GI,Number=.,Type=String,Description="Gene ID">
##INFO=<ID=LO,Number=.,Type=String,Description="Genetic location">
##INFO=<ID=FC,Number=.,Type=String,Description="Functional Consequence">
##INFO=<ID=PC,Number=.,Type=String,Description="Protein level change">
##INFO=<ID=DC,Number=.,Type=String,Description="DNA level change">
##INFO=<ID=EXON,Number=.,Type=String,Description="Exon number">
##INFO=<ID=CI,Number=.,Type=String,Description="Clinical impact">
##INFO=<ID=dbSNP,Number=.,Type=String,Description="dbSNP version">
##INFO=<ID=1KG,Number=0,Type=String,Description="Existing in the 1000 Genome phase 3 data">
##FILTER=<ID=R8,Description="IndelRepeatLength is greater than 8">
##FILTER=<ID=R8.1,Description="IndelRepeatLength of a monomer is greater than 8">
##FILTER=<ID=R8.2,Description="IndelRepeatLength of a dimer is greater than 8">
##FILTER=<ID=sb,Description="Variant strand bias high">
##FILTER=<ID=sb.s,Description="Variant strand bias significantly high (only for SNV)">
##FILTER=<ID=rs,Description="Variant with rs (dbSNP) number in a non-core gene">
##FILTER=<ID=FP,Description="Possibly false positives due to high similarity to off-target regions">
##FILTER=<ID=NC,Description="Noncoding INDELs on non-core genes">
##FILTER=<ID=lowDP,Description="low depth variant">
##FILTER=<ID=Benign,Description="Benign variant">
##contig=<ID=chrM,length=16569>
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr3,length=198295559>
##contig=<ID=chr4,length=190214555>
##contig=<ID=chr5,length=181538259>
##contig=<ID=chr6,length=170805979>
##contig=<ID=chr7,length=159345973>
##contig=<ID=chr8,length=145138636>
##contig=<ID=chr9,length=138394717>
##contig=<ID=chr10,length=133797422>
##contig=<ID=chr11,length=135086622>
##contig=<ID=chr12,length=133275309>
##contig=<ID=chr13,length=114364328>
##contig=<ID=chr14,length=107043718>
##contig=<ID=chr15,length=101991189>
##contig=<ID=chr16,length=90338345>
##contig=<ID=chr17,length=83257441>
##contig=<ID=chr18,length=80373285>
##contig=<ID=chr19,length=58617616>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr21,length=46709983>
##contig=<ID=chr22,length=50818468>
##contig=<ID=chrX,length=156040895>
##contig=<ID=chrY,length=57227415>
##contig=<ID=chrUn_GL000195v1,length=182896>
##contig=<ID=chrUn_GL000213v1,length=164239>
##contig=<ID=chrUn_GL000214v1,length=137718>
##contig=<ID=chrUn_GL000216v2,length=176608>
##contig=<ID=chrUn_GL000218v1,length=161147>
##contig=<ID=chrUn_GL000219v1,length=179198>
##contig=<ID=chrUn_GL000220v1,length=161802>
##contig=<ID=chrUn_GL000224v1,length=179693>
##contig=<ID=chrUn_GL000226v1,length=15008>
##contig=<ID=chrUn_KI270302v1,length=2274>
##contig=<ID=chrUn_KI270303v1,length=1942>
##contig=<ID=chrUn_KI270304v1,length=2165>
##contig=<ID=chrUn_KI270305v1,length=1472>
##contig=<ID=chrUn_KI270310v1,length=1201>
##contig=<ID=chrUn_KI270311v1,length=12399>
##contig=<ID=chrUn_KI270312v1,length=998>
##contig=<ID=chrUn_KI270315v1,length=2276>
##contig=<ID=chrUn_KI270316v1,length=1444>
##contig=<ID=chrUn_KI270317v1,length=37690>
##contig=<ID=chrUn_KI270320v1,length=4416>
##contig=<ID=chrUn_KI270322v1,length=21476>
##contig=<ID=chrUn_KI270329v1,length=1040>
##contig=<ID=chrUn_KI270330v1,length=1652>
##contig=<ID=chrUn_KI270333v1,length=2699>
##contig=<ID=chrUn_KI270334v1,length=1368>
##contig=<ID=chrUn_KI270335v1,length=1048>
##contig=<ID=chrUn_KI270336v1,length=1026>
##contig=<ID=chrUn_KI270337v1,length=1121>
##contig=<ID=chrUn_KI270338v1,length=1428>
##contig=<ID=chrUn_KI270340v1,length=1428>
##contig=<ID=chrUn_KI270362v1,length=3530>
##contig=<ID=chrUn_KI270363v1,length=1803>
##contig=<ID=chrUn_KI270364v1,length=2855>
##contig=<ID=chrUn_KI270366v1,length=8320>
##contig=<ID=chrUn_KI270371v1,length=2805>
##contig=<ID=chrUn_KI270372v1,length=1650>
##contig=<ID=chrUn_KI270373v1,length=1451>
##contig=<ID=chrUn_KI270374v1,length=2656>
##contig=<ID=chrUn_KI270375v1,length=2378>
##contig=<ID=chrUn_KI270376v1,length=1136>
##contig=<ID=chrUn_KI270378v1,length=1048>
##contig=<ID=chrUn_KI270379v1,length=1045>
##contig=<ID=chrUn_KI270381v1,length=1930>
##contig=<ID=chrUn_KI270382v1,length=4215>
##contig=<ID=chrUn_KI270383v1,length=1750>
##contig=<ID=chrUn_KI270384v1,length=1658>
##contig=<ID=chrUn_KI270385v1,length=990>
##contig=<ID=chrUn_KI270386v1,length=1788>
##contig=<ID=chrUn_KI270387v1,length=1537>
##contig=<ID=chrUn_KI270388v1,length=1216>
##contig=<ID=chrUn_KI270389v1,length=1298>
##contig=<ID=chrUn_KI270390v1,length=2387>
##contig=<ID=chrUn_KI270391v1,length=1484>
##contig=<ID=chrUn_KI270392v1,length=971>
##contig=<ID=chrUn_KI270393v1,length=1308>
##contig=<ID=chrUn_KI270394v1,length=970>
##contig=<ID=chrUn_KI270395v1,length=1143>
##contig=<ID=chrUn_KI270396v1,length=1880>
##contig=<ID=chrUn_KI270411v1,length=2646>
##contig=<ID=chrUn_KI270412v1,length=1179>
##contig=<ID=chrUn_KI270414v1,length=2489>
##contig=<ID=chrUn_KI270417v1,length=2043>
##contig=<ID=chrUn_KI270418v1,length=2145>
##contig=<ID=chrUn_KI270419v1,length=1029>
##contig=<ID=chrUn_KI270420v1,length=2321>
##contig=<ID=chrUn_KI270422v1,length=1445>
##contig=<ID=chrUn_KI270423v1,length=981>
##contig=<ID=chrUn_KI270424v1,length=2140>
##contig=<ID=chrUn_KI270425v1,length=1884>
##contig=<ID=chrUn_KI270429v1,length=1361>
##contig=<ID=chrUn_KI270435v1,length=92983>
##contig=<ID=chrUn_KI270438v1,length=112505>
##contig=<ID=chrUn_KI270442v1,length=392061>
##contig=<ID=chrUn_KI270448v1,length=7992>
##contig=<ID=chrUn_KI270465v1,length=1774>
##contig=<ID=chrUn_KI270466v1,length=1233>
##contig=<ID=chrUn_KI270467v1,length=3920>
##contig=<ID=chrUn_KI270468v1,length=4055>
##contig=<ID=chrUn_KI270507v1,length=5353>
##contig=<ID=chrUn_KI270508v1,length=1951>
##contig=<ID=chrUn_KI270509v1,length=2318>
##contig=<ID=chrUn_KI270510v1,length=2415>
##contig=<ID=chrUn_KI270511v1,length=8127>
##contig=<ID=chrUn_KI270512v1,length=22689>
##contig=<ID=chrUn_KI270515v1,length=6361>
##contig=<ID=chrUn_KI270516v1,length=1300>
##contig=<ID=chrUn_KI270517v1,length=3253>
##contig=<ID=chrUn_KI270518v1,length=2186>
##contig=<ID=chrUn_KI270519v1,length=138126>
##contig=<ID=chrUn_KI270521v1,length=7642>
##contig=<ID=chrUn_KI270522v1,length=5674>
##contig=<ID=chrUn_KI270528v1,length=2983>
##contig=<ID=chrUn_KI270529v1,length=1899>
##contig=<ID=chrUn_KI270530v1,length=2168>
##contig=<ID=chrUn_KI270538v1,length=91309>
##contig=<ID=chrUn_KI270539v1,length=993>
##contig=<ID=chrUn_KI270544v1,length=1202>
##contig=<ID=chrUn_KI270548v1,length=1599>
##contig=<ID=chrUn_KI270579v1,length=31033>
##contig=<ID=chrUn_KI270580v1,length=1553>
##contig=<ID=chrUn_KI270581v1,length=7046>
##contig=<ID=chrUn_KI270582v1,length=6504>
##contig=<ID=chrUn_KI270583v1,length=1400>
##contig=<ID=chrUn_KI270584v1,length=4513>
##contig=<ID=chrUn_KI270587v1,length=2969>
##contig=<ID=chrUn_KI270588v1,length=6158>
##contig=<ID=chrUn_KI270589v1,length=44474>
##contig=<ID=chrUn_KI270590v1,length=4685>
##contig=<ID=chrUn_KI270591v1,length=5796>
##contig=<ID=chrUn_KI270593v1,length=3041>
##contig=<ID=chrUn_KI270741v1,length=157432>
##contig=<ID=chrUn_KI270742v1,length=186739>
##contig=<ID=chrUn_KI270743v1,length=210658>
##contig=<ID=chrUn_KI270744v1,length=168472>
##contig=<ID=chrUn_KI270745v1,length=41891>
##contig=<ID=chrUn_KI270746v1,length=66486>
##contig=<ID=chrUn_KI270747v1,length=198735>
##contig=<ID=chrUn_KI270748v1,length=93321>
##contig=<ID=chrUn_KI270749v1,length=158759>
##contig=<ID=chrUn_KI270750v1,length=148850>
##contig=<ID=chrUn_KI270751v1,length=150742>
##contig=<ID=chrUn_KI270752v1,length=27745>
##contig=<ID=chrUn_KI270753v1,length=62944>
##contig=<ID=chrUn_KI270754v1,length=40191>
##contig=<ID=chrUn_KI270755v1,length=36723>
##contig=<ID=chrUn_KI270756v1,length=79590>
##contig=<ID=chrUn_KI270757v1,length=71251>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  TN21-116928
chr1    6197766 rs2294714   T   C   0.0 rs  DP=1931;TI=NM_000983.3;GI=RPL22;FC=Silent;DC=c.13-10A>G;LO=INTRON;EXON=2    GT:VF:AD:SB:SA:SP:SO    0/1:0.483:999,932:883/1048:451/548,432/500:2:0.0
chr1    7677739 rs3011926   A   G   0.0 rs  DP=851;TI=NM_015215.3;GI=CAMTA1;FC=Silent;DC=c.2914+6A>G;LO=INTRON;EXON=11  GT:VF:AD:SB:SA:SP:SO    1/1:0.999:1,850:309/542:0/1,309/541:0:INF
chr1    9724365 rs11121484  C   T   0.0 rs  DP=1664;TI=NM_005026.4;GI=PIK3CD;FC=Synonymous;PC=Y936=;DC=c.2808C>T;LO=EXON;EXON=22    GT:VF:AD:SB:SA:SP:SO    0/1:0.455:907,757:842/822:448/459,394/363:6:0.0
chr1    10258594    rs12402052  C   G   0.0 rs  DP=1400;TI=NM_015074.3;GI=KIF1B;FC=Synonymous;PC=A95=;DC=c.285C>G;LO=EXON;EXON=4    GT:VF:AD:SB:SA:SP:SO    0/1:0.462:751,647:678/722:375/376,302/345:6:0.1
chr1    10295776    rs3753037   C   T   0.0 rs  DP=1050;TI=NM_015074.3;GI=KIF1B;FC=Silent;DC=c.1639+10C>T;LO=INTRON;EXON=17 GT:VF:AD:SB:SA:SP:SO    0/1:0.477:549,501:395/655:211/338,184/317:2:0.0
chr1    10375266    rs11121552  C   A   0.0 rs  DP=1487;TI=NM_015074.3;GI=KIF1B;FC=Synonymous;PC=T1721=;DC=c.5163C>A;LO=EXON;EXON=46    GT:VF:AD:SB:SA:SP:SO    0/1:0.484:767,719:870/617:443/324,427/292:3:0.0
chr1    11121270    rs11121691  C   T   0.0 .   DP=1092;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=L2303=;DC=c.6909G>A;LO=EXON;EXON=49 GT:VF:AD:SB:SA:SP:SO    0/1:0.483:565,527:664/428:346/219,318/209:1:0.0
chr1    11130589    rs2275527   G   A   0.0 .   DP=1048;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=S1851=;DC=c.5553C>T;LO=EXON;EXON=39 GT:VF:AD:SB:SA:SP:SO    0/1:0.509:515,533:572/476:275/240,297/236:3:0.0
chr1    11145001    rs1057079   C   T   0.0 .   DP=1392;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=A1577=;DC=c.4731G>A;LO=EXON;EXON=33 GT:VF:AD:SB:SA:SP:SO    0/1:0.511:681,711:703/689:352/329,351/360:4:0.0
chr1    11228701    rs1064261   G   A   0.0 .   DP=1859;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=N999=;DC=c.2997C>T;LO=EXON;EXON=19  GT:VF:AD:SB:SA:SP:SO    0/1:0.490:948,911:1059/800:536/412,523/388:2:0.0
chr1    11241657    rs1135172   A   G   0.0 .   DP=1643;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=D479=;DC=c.1437T>C;LO=EXON;EXON=10  GT:VF:AD:SB:SA:SP:SO    0/1:0.512:801,842:793/850:376/425,417/425:5:0.0

Unfortunately, when I try to utilize the following Python code to validate my VCF file:
import pysam

def validate_vcf(vcf_file_path):
    try:
        # Open the VCF file
        vcf_file = pysam.VariantFile(vcf_file_path)

        # Iterate through the records to check if they can be read
        for record in vcf_file:
            pass

        print("VCF file is valid.")
    except Exception as e:
        print(f"Error validating VCF file: {str(e)}")

# Path to the VCF file you want to validate
vcf_file_path = "path/to/yourfile.vcf"

validate_vcf(vcf_file_path)

I receive the following error:
[W::bcf_hrec_check] Invalid tag name: "1KG"
[E::vcf_parse_format] Invalid character '/' in 'SA' FORMAT field at chr1:6197766

Does anyone know what I can do to fix this?

It looks like the validator is working, in that it's producing messages about validation. Why do you need this validator? What are you trying to use the VCF file for?

Given issues I've had in various programming languages that use numbers at the start of a variable, the tag issue might be related to it starting with a numeric value. Is it possible to change that from 1KG to something like ph3_1KG?

According to the information in the VCF header, the SA field should have 4 values:
##FORMAT=<ID=SA,Number=4,Type=Integer,Description="Number of 1) forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles, used in variant calling">

While that location does have four numerical values in the actual data for that variant location, it seems like the validation code isn't expecting the forward and reverse alleles to be separated by /:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  TN21-116928
chr1    6197766 rs2294714   T   C   0.0 rs  DP=1931;TI=NM_000983.3;GI=RPL22;FC=Silent;DC=c.13-10A>G;LO=INTRON;EXON=2    GT:VF:AD:SB:SA:SP:SO    0/1:0.483:999,932:883/1048:451/548,432/500:2:0.0
                                                                                ^^                                     ^^^^^^^^^^^^^^^

Whether or not these errors are actually important will depend on the program that this VCF file is used for.


text:
WDL dependencies
I thought that when you linked tasks in a WDL workflow e.g.
call task1 {
    input: in = file1
}

call task2 {
    input: in = task1.out
}

That the first time you run it task1 and task2 would run, but the second time you run it nothing would run because the inputs haven't changed i.e. like GNU make
Also, if you delete the output of task2 and run it then only task2 would run as the output from task1 is still there and hasn't changed.
Basically, I am trying to use WDL like GNU make to only re-run the parts of my workflow that need rerunning depending on which inputs have changed.
I suspect I have misunderstood and WDL cannot be used like this?

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

Welcome to the site and thanks for the question. Could you provide a reproducible example for the question? This does help understand what is sought. By the way, WDL is not hugely popular on this site, however NextFlow however is popular.

This is dependent on the execution engine you're using. If you are using cromwell or miniwdl then you can try looking up "call caching" in their documentation. Both of those execution engines should support it.


text:
Issue creating CNV plot from WES data
I received Whole Exome Sequencing data from an NGS company (CARIS, specifically). I received R1 and R2 FASTQ files, a BAM file aligned to hg38, and a VCF file.
I used CNVPytor to create a CNV plot (Manhattan plot), by following this example Jupyter Notebook here: https://github.com/abyzovlab/CNVpytor/blob/master/examples/PythonLibraryGuide.ipynb
However, when I run this code on my data, I get the following graph:

Furthermore, I tried to use CNVKit and used the following steps/code:
cnvkit.py batch TN21-116928.DNA.bam --normal -m hybrid --fasta hg38.fa --targets targets.bed --output-reference my_reference.cnn

Where "targets.bed" was a file of the following form:
chr1    33306766    33321098    A3GALT2
chr22   42692121    42721298    A4GALT
chr3    138123713   138132390   A4GNT
chr12   53307456    53324864    AAAS
chr12   125065434   125143333   AACS
chr3    151814073   151828488   AADAC

The graph created from this is the following:

Does anyone know where I am going wrong? Any pointers?

cnvkit.py scatter -s Sample.cn{s,r} -c chr7 maybe? Anyway its a cool question (except for the results)


text:
What are the fastest way of feature selection in R for large dataset in Machine Learning?
I have about 100 data sets with over 100,000 rows and 27 predictor variables and multiple response variable on each dataset. I want to select the most essential Features for these datasets. I used Boruta() from the package Boruta, but it takes an extremely long time (about 3 hours for a single sample). Could someone kindly tell me the quickest technique to select key features from a large dataset?
This is the extension from my previous questions:

https://stackoverflow.com/questions/76824529/how-to-do-regression-modelling-in-r-multiple-one-at-a-time-to-multiple-column
How to get coefficient tables from multiple regression model result!

Your second question has nothing to do with machine learning, it's a simple R question - you ran into problem because you did not use the apply family of functions properly.

Right, I thought that too. Really an analyst uses either ML or GLM - not both. Within ML it's common to compare different methods but not move outside ML (unless they want to do deep learning). Its possible "linear regression" ML might be intended. ML is very, very different to GLM, i.e. the method and system to test model.

Feature selection: the fastest thing moving (i.e. the simplest) is naive Bayes ML. I don't use it but KNN (K-nearest neighbor) I do use, which is also computationally simple, i.e. fast. I think it gives better results, i.e. the best results for a simple algorithm.
R ML has improved decision trees BTW - quite popular, but you'd never use this in Python.


text:
The current status of human genome annotation
I was wondering if you can suggest something on the current status of the human genome annotation efforts. I know about ENSEMBL/HAVANA/GENCODE, and I've read the relevant papers; however, I was wondering if RefSeq or UCSC are still actively developed, and are there any competing efforts on T2T assembly, human pan-genome, etc etc?
I know there are quite a few projects that generate both short and long read RNA-seq and other "evidence" from lots of cell lines etc. I guess I was looking for a good and up-to-date perspective or a review type of article.
Thank you in advance!

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Welcome to the site @AlexanderPredeus. Thanks for your question, there are a number of different questions here and the format we request is focus on a single issue as the basis of the question. You can always ask multiple questions. What we want to avoid is a discussion forum. Please do consider editing your question.


text:
IgBLAST Error: Germline annotation database could not be found in [internal_data] directory
I'm in the following situation and I came here to see if I could get some help.
I'd like to use IgBLAST for an analysis on some Alpaca antibody data, so I received some FASTA sequences from a colleague, imgt_gapped_alpaca_IGHV.fasta, imgt_gapped_alpaca_IGHD.fasta, etc.
I used the edit_imgt_file.pl tool to appropriately structure each .fasta file, and to make a BLAST database for each (V, D, J) segment:
# Process V segment sequences with edit_imgt_file
edit_imgt_file.pl \
~/VHH/data/imgt-data/imgt_gapped_alpaca_IGHV.fasta > ~/ncbi-igblast-1.21.0/internal_data/alpaca/imgt_alpaca_ig_v.fasta

makeblastdb -parse_seqids -dbtype nucl \
    -in ~/ncbi-igblast-1.21.0/internal_data/alpaca/imgt_alpaca_ig_v.fasta \
    -out ~/ncbi-igblast-1.21.0/internal_data/alpaca/alpaca_V

After doing that for each of the segments, I went to put all the pieces together and perform the following query:
igblastn \
    -germline_db_V ~/ncbi-igblast-1.21.0/internal_data/alpaca/alpaca_V \
    -germline_db_D ~/ncbi-igblast-1.21.0/internal_data/alpaca/alpaca_D \
    -germline_db_J ~/ncbi-igblast-1.21.0/internal_data/alpaca/alpaca_J \
    -auxiliary_data ~/ncbi-igblast-1.21.0/optional_file/alpaca_gl.aux \
    -domain_system imgt -ig_seqtype Ig -organism alpaca \
    -outfmt '7 std qseq sseq btop' \
    -query ~/VHH/data/ALPACA1-BC1.fastq \
    -out ~/VHH/data/igblast-out/ALPACA1-BC1.fmt7

Unfortunately, the following error is returned to me:
BLAST query/options error: Germline annotation database alpaca/alpaca_V could not be found in [internal_data] directory
Please refer to the BLAST+ user manual.

I've searched the internet high and low, but this appears to be a bit of a niche problem, and the BLAST+ user manual promised a little more than it was able to deliver in terms of helpful suggestions.
Any questions/comments/ideas/diagnostics would be very welcome-- if you see any low hanging fruit, please let me know-- thank you so much for your time!

What you have looks like the right idea to me.  (I use igblastn regularly for human and rhesus macaque antibodies with custom DB+aux files, but not with a custom organism label itself.)  What files do you see in your alpaca db dir for V/D/J?  Are both igblastn and makeblastdb definitely coming from igblast 1.21.0?  (Any chance something's mismatched between those two binaries, I wonder?)  And where did your alpaca auxiliary file come from?  You could maybe drop the aux part to simplify things while troubleshooting.

Actually I get the same error if I just change my human or rhesus commands to some arbitrary organism name.  Did you also create the FWR/CDR file like it talks about here?  I guess all of my questions really boil down to, what's in your ~/ncbi-igblast-1.21.0/internal_data/alpaca directory?


text:
enrichment of DE genes within wide regions not pathways, hypergeometric test or GSEA?
I haven't found any relevant posts, so here is my question (it's more a question on statistics I believe):
I have the deseq2 output of DE genes( around 25000 genes) , and among them 1000 are significantly differentiated(and I want to check them). Now, I group these 25000 genes by the epigenetic states of my data(i divide the full genome into 3 groups and I also have 4th group which are regions whose group I couldn’t identify due to lack of coverage etc ). So now I have my 25,000 genes where roughly 12000 belong to group a, 8000 belong to group 4 unidentified, 4500 belong to group b and 500 belong to group c. Now, what i do , from these 25000 genes I define my 4 gene sets where each group is one term so 4 in total, then i take my 1000 DEGs , i rank them and then I use prerank approach to get the enrichment results. I also hesitate to just stick ot hypergeometric test instead, but then not sure what I should consider in this case.
My question is: How can I check if there are more differentially expressed genes in one group compared to others? I am particularly interested in group A, i expect it to have more DE genes (i don't know if upregulated or downregulated or both mixed).
So far I've decided to go with hypergeometric test, but if my hypothesis is that one group is overrepresented with DE than others, i'm not sure it's a good approach; I tried GSEA too but because i have very few large gene sets ( where i define by a gene set all the genes that fall into one group), i hesitate if it's a suitable approach.
I would greatly appreciate any help and tips answers and comments.  Thanks a lot everyone!!
best regards,
alexa

@M__ I've corrected it. hope it's fine now. thank you

Thanks @drew, it's still not absolutely clear but is clear enough.


text:
Get SRR numbers from a GSE number
I'm trying to programmatically download files from SRA (with prefetch in sra-tools) starting from a GSE number, using Python.
There is an explanation for using entrez-direct on the CLI, using esearch and the efetch. But I want to use Python. This is what I have now:
from Bio import Entrez
Entrez.email = 'my_name@email.me'
handle = Entrez.esearch(db='gds', query='GSE134195', term="biopython")
record = Entrez.read(handle)

which gives:
record = Entrez.read(handle)
record

{'Count': '17', 'RetMax': '17', 'RetStart': '0', 'IdList': ['200066260', '200045602', '200037064', '200030751', '100013939', '300762998', '300762997', '300762996', '300762995', '300762994', '300762993', '300762992', '300762991', '300762990', '300762989', '300762988', '300762987'], 'TranslationSet': [], 'TranslationStack': [{'Term': 'biopython[All Fields]', 'Field': 'All Fields', 'Count': '17', 'Explode': 'N'}, 'GROUP'], 'QueryTranslation': 'biopython[All Fields]'}

But what are these numbers? IdList looks useful, but there are 17 numbers where I expect 14. How do I get to SRR (SRA) files? The files I want can be found manually here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134195, click on SRA run selector (at the bottom), this bring you here: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA554267, there I can download the Accession list that I want. But I'd like to automate that process.
Thank you in advance.

I think thats easy, if I get chance I'll put together a bit of code.

That would be very nice of you, thank you.


text:
imgtrefseq vs imgtblast
In IMGT ftp-download (flat files) page, there are 2 optional fasta files (imgtblast and imgtrefseq) with their corresponding blast databases to download. Is the only difference the syntax and additional annotation in the refseq files?
Thank you!


text:
What alternatives are there to Bracken?
Bracken is the recommended post-processer of Kraken2 outputs and specifically is expected to improve relative abundance estimates. However, Bracken's Github has over 100 outstanding issues and the code seems to rely on Python 2, which is outdated. While many of the github issues are likely not due to bugs, it does seem like the software does not seem to be supported strongly at this time. I am curious what other alternatives exist. Are folks moving from kraken2/bracken to another solution for quick taxonomic hypothesis generation of metagenomics data?

If you want guaranteed support, then pay for a commercial solution ;)

Happy to consider that. What commercial solutions do you prefer or pay for?

No, Braken is the tour de force of first stage post-hoc Kraken2 analysis.
Whether the code is supported or not it corrects bias in Kraken2 population frequency distributions. I'm not aware of any theoretical concerns.
Generally, Python2 isn't a good idea, but it is overshadowed with the bigger concern of population bias, which is what Braken is there to adjust.


text:
Try to understand better how to set out an effective sample size within the MetagenomSeq package in microbiome
Cross-posted on biostars: https://www.biostars.org/p/9572744/
I am trying to understand better what number to be used and how to put it in the function MRcoefs within the Metagenomeseq package in R for calculating differentially expressed OTUs between two condition:
my code looks like:
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = 5.5)

Here I am setting the eff argument to 5.5 because I have calculated the effectiv_size using this code:
Effective_size(fit)

which result in table of effective size for each feature (OTU), and then I got the mean of these
mean(Effective_size(fit))

which was 5.5
When I run this, it produced an error:
prop outside [0,1], I understand it should be between 0-1, but why ? is not that number in eff is the mean that I calculated which bring 5.5 ?
Also when I tried to run it with .5 as the document suggested as below :
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = .5)

it run nicely, no error , but it produced the same OTUs as when I did not include eff, I was expecting that some OTUs will be removed as they were below this threshold
Could you please help me ?

Hi @MohamedSamir, firstly please put your code into coding blocks. Secondly, the package you are using is Metagenomeseq but the tag is 16S. These are two different approaches, could you please clarify which you are using?

Thanks for your advice. I mean Metagenomeseq which is used for analysing 16S microbiome data. I mean MetagenomeSeq package

Set numberEff = True
MRcoefs(fit, coef = 2, group = 3, number = 143, eff = 5.5, numberEff = True)

I suspect you've got the quantile value set, hence it will throw an error above 1.


text:
A bash script for running on a bunch of bam files
I have some bam files in this directory
/data/Continuum/WES/results/

I want to run GATK mutation calling over bam files
I googled and I realised for one function, I can do this
cd /data/Continuum/WES/vcf/

for file in *.bam ; do ./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta done;

I have five functions to run but I do not know how to include all in the script
These are my functions
./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

./gatk GetPileupSummaries -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -V af-only-gnomad.hg38.vcf.gz -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

./gatk CalculateContamination -I NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.contamination.table

./gatk GetSampleName -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name

./gatk Mutect2 -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.vcf -tumor NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name --af-of-alleles-not-in-resource 2.5e-06 --germline-resource af-only-gnomad.hg38.vcf.gz -pon 1000g_pon.hg38.vcf.gz

I thought to use | to include five functions but I do not know how to change input and output for each iteration
I tried this for the first function,
    FILES="/data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*"
for f in $FILES
do
  ./gatk CollectSequencingArtifactMetrics -I $f -O $f --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta
done

but returned this error
SAMException: Cannot read non-existent file: file:///data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*

Anyone, please help me in develop this?

There are a couple ways to approach this, one is to create output files based on the BAM input file name - the other approach is to make a new folder (based on the input BAM filename containing result files with the same name for each input. I generally use the first method, but either work.
Here is a script that may run the 5 GATK functions. It calculates a prefix to be used for derivative output files BASE, and should run your pipeline on all the BAM files within the INPUT_FOLDER.
#!/bin/bash

set -e

INPUT_FOLDER="/data/Continuum/WES/results/"

for BAM_FILE in "$INPUT_FOLDER/"*.bam
do
        BASE=$(basename "$BAM_FILE" .bam) # NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup

        # function 1
        ./gatk CollectSequencingArtifactMetrics \
                -I "$BAM_FILE" \
                -O "$BASE" \
                --FILE_EXTENSION .txt \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 2
        ./gatk GetPileupSummaries \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.table" \
                -V af-only-gnomad.hg38.vcf.gz \
                -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 3
        ./gatk CalculateContamination \
                -I "$BASE.targeted_sequencing.table" \
                -O "$BASE.targeted_sequencing.contamination.table" \

        # function 4
        ./gatk GetSampleName \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.sample_name"

        # function 5
        ./gatk Mutect2 \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
                -I "$BAM_FILE" \
                -O "$BASE".vcf \
                -tumor "$BASE.targeted_sequencing.sample_name" \
                --af-of-alleles-not-in-resource 2.5e-06 \
                --germline-resource af-only-gnomad.hg38.vcf.gz \
                -pon 1000g_pon.hg38.vcf.gz
done

It may need some adjustments - I did not test this, and used the GATK incantations you provided (thank you for including all of your steps!)
It includes the addition of the command set -e, which depending on GATKs return status may stop the script if any command fails - otherwise the script will happily run through a failed command.
You also may need to adjust the location of your output, and the directory the script is in upon execution, as you are running GATK from the current directory (./gatk)
In order to make these long lines readable, I use the method of placing a \ at the end of a line - this allows the command to wrap to the next line. Be sure that the \ character is at the end of your line for this to function properly.

Wow, cool, thanks.

OMG Thanks a million the script is running

Sorry @Scot, as the Mutect2 part of the script takes 15 hours in average for each sample, does this mean that for function 5, script will wait for the first sample to being completed by Mutect2 then goes for the second sample or Mutect2 will run of all of 18 samples in parallel? In the first case, the script will take 18 days?

Shell scripts are run line by line. If a line takes 15 hours, the interpreter will not move to the next line before that line is done. Loops are merely instructions to run a set of lines multiple times in series, so there's no parallel action happening here. You should ideally write a snakemake script with a rule for each GATK step then run a bunch of snakemake processes in parallel using whatever parallelization means are available to you.

Hi @Angel - yes, the current version of this script runs sequentially, and is not parallelized. If function5 Mutect2 takes a long time, then it may be worthwhile to pull it out of the loop, and run it upon completion of all other steps, and addition, explore if it can be run in parallel.

You really don't want to use CAPS for shell variable names. Since, by convention, global environment variables are capitalized, it is bad practice to also use caps for your own, local variables because this can lead to naming collisions and hard to find bugs.


text:
MergeBamAlignment error
I doing the alignment of samples following the GATK pipeline, and doing the MergeBamAlignment,like this:
MergeBamAlignment \
  -ALIGNED $path/file.unsorted.bam \
  -UNMAPPED $path/file.unmapped.bam \
  -O $path/file.merged.bam \
  -R $references/GRCh38.primary_assembly.genome.fa \
  -SO coordinate \
  -TMP_DIR $path/tmp/

I got this error:
Exception in thread "main" java.lang.IllegalArgumentException: Do not use this function to merge dictionaries with different sequences in them. Sequences must be in the same order as well. Found [] and [chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr20, chr21, chr22, chrX, chrY, chrM, GL000008.2, GL000009.2, GL000194.1, GL000195.1, GL000205.2, GL000208.1, GL000213.1, GL000214.1, GL000216.2, GL000218.1, GL000219.1, GL000220.1, GL000221.1, GL000224.1, GL000225.1, GL000226.1, KI270302.1, KI270303.1, KI270304.1, KI270305.1, KI270310.1, KI270311.1, KI270312.1, KI270315.1, KI270316.1, KI270317.1, KI270320.1, KI270322.1, KI270329.1, KI270330.1, KI270333.1, KI270334.1, KI270335.1, KI270336.1, KI270337.1, KI270338.1, KI270340.1, KI270362.1, KI270363.1, KI270364.1, KI270366.1, KI270371.1, KI270372.1, KI270373.1, KI270374.1, KI270375.1, KI270376.1, KI270378.1, KI270379.1, KI270381.1, KI270382.1, KI270383.1, KI270384.1, KI270385.1, KI270386.1, KI270387.1, KI270388.1, KI270389.1, KI270390.1, KI270391.1, KI270392.1, KI270393.1, KI270394.1, KI270395.1, KI270396.1, KI270411.1, KI270412.1, KI270414.1, KI270417.1, KI270418.1, KI270419.1, KI270420.1, KI270422.1, KI270423.1, KI270424.1, KI270425.1, KI270429.1, KI270435.1, KI270438.1, KI270442.1, KI270448.1, KI270465.1, KI270466.1, KI270467.1, KI270468.1, KI270507.1, KI270508.1, KI270509.1, KI270510.1, KI270511.1, KI270512.1, KI270515.1, KI270516.1, KI270517.1, KI270518.1, KI270519.1, KI270521.1, KI270522.1, KI270528.1, KI270529.1, KI270530.1, KI270538.1, KI270539.1, KI270544.1, KI270548.1, KI270579.1, KI270580.1, KI270581.1, KI270582.1, KI270583.1, KI270584.1, KI270587.1, KI270588.1, KI270589.1, KI270590.1, KI270591.1, KI270593.1, KI270706.1, KI270707.1, KI270708.1, KI270709.1, KI270710.1, KI270711.1, KI270712.1, KI270713.1, KI270714.1, KI270715.1, KI270716.1, KI270717.1, KI270718.1, KI270719.1, KI270720.1, KI270721.1, KI270722.1, KI270723.1, KI270724.1, KI270725.1, KI270726.1, KI270727.1, KI270728.1, KI270729.1, KI270730.1, KI270731.1, KI270732.1, KI270733.1, KI270734.1, KI270735.1, KI270736.1, KI270737.1, KI270738.1, KI270739.1, KI270740.1, KI270741.1, KI270742.1, KI270743.1, KI270744.1, KI270745.1, KI270746.1, KI270747.1, KI270748.1, KI270749.1, KI270750.1, KI270751.1, KI270752.1, KI270753.1, KI270754.1, KI270755.1, KI270756.1, KI270757.1].
    at htsjdk.samtools.SAMSequenceDictionary.mergeDictionaries(SAMSequenceDictionary.java:340)
    at picard.sam.SamAlignmentMerger.getDictionaryForMergedBam(SamAlignmentMerger.java:197)
    at picard.sam.AbstractAlignmentMerger.mergeAlignment(AbstractAlignmentMerger.java:384)
    at picard.sam.SamAlignmentMerger.mergeAlignment(SamAlignmentMerger.java:181)
    at picard.sam.MergeBamAlignment.doWork(MergeBamAlignment.java:368)
    at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:309)
    at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
    at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)

srun: error: compute-12: task 0: Exited with exit code 1

Before I also confirmed that the files that I'm trying to merge, are with the same reference sequences. So I don't understand what could be the problem.

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

Welcome @RitaSoares. Firstly is this a human genome? Secondly, what have you done unto this point, ie. how did you create the different alignments, especially what were being used?

Not sure if this is related, but do you get the same error after sorting your aligned bam file? Also, just checking that you aligned the reads to the same reference ($references/GRCh38.primary_assembly.genome.fa)

@jared_mamrot yes but i got a different error after sorting the aligned bam files, with SortSamSpark: "org.broadinstitute.hellbender.tools.spark.pipelines. SortSamSpark done. Elapsed time: 0.07 minutes.Runtime.totalMemory()=2854748160
htsjdk.samtools.util.RuntimeEOFException: Premature EOF. Expected 4 but only received 0; BinaryCodec in readmode; streamed file (filename not available)". But according to this error that happend to me before, normally i put in the code of MergeBamAlignment the argument "SO coordinate"and  this corrects my problem but with this samples doesn't happen that.

Sounds to me like your sample might be corrupted/truncated. What does picard ValidateSamFile output for that sample?

I also think that is the case... i will try to do ValidateSamFile to make sure. Thank for your answer


text:
Feature-selection microbiome data to build a model to predict cancer
I want to build a machine learning model to predict colorectal cancer based on 16S rRNA microbiome data (stool samples).
I have filtered the data using filtering approach by Duvallet(removing samples with fewer than 100 reads and OTUs with fewer than 10 reads,as well as OTUs which were present in fewer than 1% within the study) What do I have to do after the preprocessing step? How can I do the feature selection for the phenotype data and the microbiome data?

Welcome. Thanks for your question. I'm afraid this isn't how the site works. What we answer is specific technical stuff rather than advice on a project design. Preprocessing -> output of ML is an entire project. My advice is to focus on a single question supply and as much information as you can. For example is the microbiome the intestinal trait, is this bowel cancer? How are you sampling the microbiome, 16S, whole genome? You say OTUs - is that optimal taxononic units - this infers you've produced a tree.

I think the question needs improving. As the question stands pre-processing to feature selection ML (that's whats being asked).

put the species diversity into a pandas dataframe the last column is the target (cancer vs healthy).
assess a variety of ML algorithms (because they give different answers). Naive bayes, linear regression and K-nearest neighbour (KNN) are simple algorithms, easy to run and at this level fine.
each have a "feature selection" - all features will be positive.

At a personal level I like KNN, but avoid LR and NB.
I personally would simply suggest unsupervised learning.


text:
How to loop over certain files within a directory that belongs to certain samples
This question was also asked on Biostars
I have two samples where the sequencing was done on 3-lanes paired end reads. The two samples are A44943_1 and A44944_1.
I would like to generate a for in loop to merge the R1 files of sample A44943_1 from the 3-lanes,  and then do the same for R2 files of this sample, then I want to merge R1 files of sample A44944_1 from the 3-lanes and then do the same for R2 files of this sample. So I thought about for in loop. My fastq files looks like:
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L001_R1_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L001_R2_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L002_R1_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L002_R2_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H7YKTDRXY_CGCAACTA-GAATCCGA_L003_R1_trimmed.fastq
R0480-S0001_182602S1P1_A44943_1_H7YKTDRXY_CGCAACTA-GAATCCGA_L003_R2_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L001_R1_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L001_R2_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L002_R1_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L002_R2_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H7YKTDRXY_CACAGACT-TGGTACAG_L003_R1_trimmed.fastq
R0480-S0002_182602S1P2_A44944_1_H7YKTDRXY_CACAGACT-TGGTACAG_L003_R2_trimmed.fastq

Could you advise me on  a for in loop to do this task?
Thanks

You can certainly use a for loop for this if desired, but an easier method is to use the cat command - something like this:
cat R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L00[123]_R1_trimmed.fastq > R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_R1_trimmed.fastq
cat R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_L00[123]_R2_trimmed.fastq > R0480-S0001_182602S1P1_A44943_1_H77YFDRXY_CGCAACTA-GAATCCGA_R2_trimmed.fastq
cat R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L00[123]_R1_trimmed.fastq > R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_R1_trimmed.fastq
cat R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_L00[123]_R2_trimmed.fastq > R0480-S0002_182602S1P2_A44944_1_H77YFDRXY_CACAGACT-TGGTACAG_R2_trimmed.fastq

Holy s**t, have I been missing out all along? TIL, thank you @Scot!

I can relate, @RamRS!  I love to learn a new trick... and bash is full of them. :)

This is news to me too so I had to look this one up.  Evidently it's part of bash's "pathname expansion" feature set, just like * (and so, distinct from "brace expansion" in that it depends on what files are present; so like, cat [123].txt expands to cat 1.txt 2.txt 3.txt if those are present, but is just cat [123].txt if not, versus {1..3}.txt which always expands).  Very handy, thanks Scot.