text:
MEGA11 gives the error "stop codons are found in the translated sequences", even though there are not!
I am using MEGA11 MEGA11: Molecular Evolutionary Genetics Analysis Version 11    ,   Mol. Biol. Evol. 38(7):3022–3027   to align protein coding sequences that I downloaded from KEGG database. I am doing MUSCLE alignment and I choose the "align codons" option. The problem that occurs has to do with the error "stop codons are found in the translated sequences" that pops as an warning during the process. When I check the sequences I cannot find any stop codons except the one at the end, as expect. I don't know how to interpret the message. Does anyone has any idea what is going wrong? Thank you in advance!

I think this is caused by the final stop codon. Have you tried removing the stop codon at the end?

https://www.genome.jp/kegg/pathway.html   there is genes https://www.genome.jp/kegg/genes.html

Thanks @pippo1980 I did notice this. Wasn't sure how the database could be cleanly interrogated easily however.

What is the format of your input file? Can you show us an example?

The solution is very simply, MegaX will not allow any stop codons within the "translated" sequence. The reason is simple

MEGAX's function is phylogenetic models with nucleotide or protein data
All phylogenetic protein calculations require a 20x20 matrix of amino acid frequencies, e.g. BLOSUM, Dayhoff, mtREV24 etc ... No matrix has ever calculated stop codon frequency (well NCBI have and sort of obviously always '1' or transformed '1' [the definition of a gene is it has one stop codeon])

Thus MEGAX will not permit a single stop codon, even the terminal codon of the gene, because it can't calculate phylogenetic information from them (doesn't fit the model). Thus the "getter" to pick up bugs occurs right at that beginning.
Two solutions:

Don't translate the nucleotide sequence on import - just stick with nucleotides (there are no codon models in MEGAX).
Remove the terminal 3 resides

     perl -p -i -e 's/(^[AGCT-NRY]{3}.*)[AGCT-NRY]{3}$/$1/g' alignment.fa

Note the above code is not tested and intended only for fasta sequence or MEGA input format - don't use it on Clustal or Phylip it will be disastrous, i.e. interleaved format. IN ADDITION make sure there are no carriage returns within the contiguous sequence (again the result is disastrous). Finally, note I forget the {3} stuff (its ages since I wrote perl) and should include all the degenerate nucleotide code beyond NRY, e.g. MB etc ...
Really trimming the 3' last 3 nucleotides is a good use of Biopython - no risk of errors regardless of alignment. Trimal can probably do it BTW or seqkit.

I'm still interested in how you interrogated KEGG gene DB BTW

If you just want to remove the last three character from every line (but why are we assuming single-line sequences here, and shouldn't we skip the fast headers?), you can use perl -i -pe 's/...$/' file. Your command will remove the last three from every line, irrespective of whether they're actually the last three of the gene.

If this is fasta then the header > or # would also have 3 characters removed. Using ^[AGCT-NRY]{3} ensures only the sequence is altered. However, many formats of MEGA and fasta using a flxed column width of say 80 characters per line. Personally, I'd use BioPython within their sequence object - which is brilliant.

I think that is it is possible to just leave it like the way it is and just choose the "ignore" option that pops. When I do that the programm goes on, however I don't know if there are any consequences for my alignment. By the way thank you everyone for your answers! I do not think that it is a big deal however I wasn't sure what caused the warning so I wanted your opinion about it. Thanks again!


text:
Aligning sequences with multiple genetic codes!
I am doing a project on duplicated genes and I have a major difficulty on how to align sequences that use different genetic codes. I work with fasta files that contain sequences of protein coding genes, every fasta file includes genes from many species that are orthologs/paralogs. These genes even in the same species may use different genetic codes. My goal is to align the sequences, codon alignment [the process of codon alignment considers the reading frame of the translated protein and adjusts the nucleotide alignment so that the protein alignment stays in frame.] , and afterwards to create N-J trees. For the whole process I work on MEGA11. I cannot organize the species-genes into groups with common genetic codes because I need all of them in the tree. One option, to my mind at least, would be the alignment not in codon but in DNA,or another option is to align their protein sequences or lastly a third option could be the change of the sequences in order to get rid of the stop codons and afterwards to use one common code for everything, by doing that I know that I will lose some info but at least I can work.Do you have any idea how to overcome this obstacle? Do you agree with any of this options? Thank you in advance!

Can you say a little more about how diverse the organisms that you're studying are? Hominins? Carnivores? Mammals? Vertebrates? Eumetazoa? Eukarya? Scale matters! the more distantly related your organisms are (including divergent genetic codes), the more you would want to use a more stable data type such as protein sequence.

Why is the genetic code relevant? There are not that many genetic codes out there, and most differences are around stop codons as you point out, but if you are aligning nucleotides, it shouldn't really matter. What does "not in codon but in DNA" mean? And yes, aligning proteins would make sense, did you try that?

I though the main differences were the relative abundancies of different codon triplets for the same aa i think they call it codon usage

@pippo1980 no, that's a different thing. There are a few different genetic codes used by some species, and even our own mitochondria. See https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG2.

Cool , how they were able to identify them ? From genetic analysis or native protein sequencing , maybe by mass spec ? Or trna and synthases study ??

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8629427/#:~:text=The%20first%20alternative%20genetic%20codes,obtained%20by%20direct%20protein%20sequencing. A computational screen for alternative genetic codes in over 250,000 genomes , here I got some background info, thanks @terdon

Its probably easy but want to avoid an explaination because as other site members point out the question lacks key details and there's risk of a "thats not what I meant" response.

@MaximilianPress I am working with every species listed on the Plants, Protists and Fungi subcategories mentioned here https://www.genome.jp/kegg/catalog/org_list.html ,these are the organims that KEGG database has.

@pippo1980  To give a better picture of my project I am working on paralog genes as they are mentioned on KEGG, there is a specific option for that on KEGGs website. After I download the sequences I align them on MEGA11 and I create N-J trees. Every tree constits of genes (paralogs+orthologs) that belong on the same KO (ortholog code). My rpoblem is that many of these paralog genes are found on different positions (mitoch. or nucleus) so I cannot align them as codons because in most cases they use differerent genetic code. My question is how can I overcome this difficulty.

@terdon I am working with gene coding sequences, do you think that if I align them as non-coding sequences there would be a major diffenence in my alignment? By "not in codon but in DNA" I mean what exactly I just wrote down here, I mentioned it like that because this is how it is mentioned on the MEGA11 program. Here I give you an example of 2 paralog genes that are located on differnent positions (mitoch. and nucleus), this in not a rare ocassion for my dataset. https://www.genome.jp/entry/taes:123096964 https://www.genome.jp/entry/taes:34688768

My sense is that you are overthinking it- just use protein sequence!

Any codon alignment tool that remove this genetic code differences , keep in mind I dont know about how much different they could be, like sterilizes the codon differences removing them from the scores or treat them as identical ? Not sure this make any sense but maybe ....

I am having a very similar issue, however mine are multiple sequences of the same specie. I have a mixture of DNA and RNA sequences so although it appears to aligning fine the proteins produced are not correct. I’m at a loss on how to fix this so if you managed to figure out what to do please let me know.

source: Stack Exchange

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

sort: votes


text:
Discrepancy in Read Counts Between FastQ and BAM Files in Adapter-Trimmed Pipeline
In a FastQ to BAM pipeline where only adapter trimming is performed, I've noticed a potential discrepancy in read counts between the initial FastQ files and their resulting BAM file. Specifically, I'm seeking clarification on whether the following statement holds true:
"Total number of reads in R1 and R2 FastQ files = mapped and unmapped reads in the corresponding BAM file."
I understand that adapter trimming is the sole preprocessing step in this pipeline. However, I'm uncertain about the impact of factors such as PCR duplicates, mapping efficiency, and any potential loss during processing. The BAM file does contain PCR duplicates, and I'm wondering if these factors contribute to the observed differences in read counts.
Can someone shed light on the expected relationship between the read counts in the FastQ files and the BAM file under these conditions? Any insights into how PCR duplicates may affect the discrepancy would be particularly helpful. Additionally, if there are common reasons for the observed differences in read counts between FastQ and BAM files, I'd appreciate guidance on identifying and addressing them.
Additional information:
The pipeline used:
pre-trimming: fastqc R1.fastq.gz R2.fastq.gz -o pre_trimming
adapter trimming: fastp -i R1.fastq.gz -o trimmed_R1.fastq.gz -I R2.fastq.gz -O trimmed_R2.fastq.gz -h file.html -j file.json
alignment: bwa mem -M -t 4 -R “@RG\tID:FE_S211\tPL:ILLUMINA\tSM:A4” hg19.fa trimmed_R1.fastq.gz trimmed_R2.fastq.gz | samtools sort -o output.bam -O BAM -@ 4 -

My main question is whether you refer to the fastq files after adapter trimming or before. It seems you mean the original fastq files. In that case, any read filtering done by fastp will remove low quality or totally-adapter reads.
Note from the fastp docs that length trimming is enabled by default, so trimmed reads will be removed if they are shorter than the length setting, even if they are otherwise high quality. And if they are low quality then they will be filtered out for sure. I refer you to the docs to see what all the various default settings are. Possibly you could turn them all off and still see the initial number of reads, but I very much doubt that they would be useful for you.
UPDATE: see also @Yossi Farjoun's excellent point about alignments vs. reads- not at all the same thing! though you can learn the exact number of reads from a BAM (using e.g. samtools stats), the BAM is in fact a record of alignments. So it's important to be looking at the right number when you're doing your comparison!

"Total number of reads in R1 and R2 FastQ files = mapped and unmapped reads in the corresponding BAM file."
there might be confusion regarding the word "read".
"record" refers to an entry in the BAM file (a row in SAM)
"template" refers to the original sequence from the fastq/sequencer.
after mapping a single template may be the source of two different records due to secondary and supplementary alignments.
So the general answer to your question is: Not if you count "reads" as rows in the SAM file Yes, if you count carefully and make sure to not double-count records that have the same "query name" (the name from the fastq file)
For more information about SAM/BAM/CRAM the spec can provide hours of light reading... :-)

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


text:
How to download Sequencing data on Windows using SRA toolkit?
I have downloaded and installed the SRA toolkit, but there seems to be no online material on how to download SRR files using the SRA toolkit on Windows.
All online materials are Ubuntu or Linux based.
Is there some obvious resource that I am missing?

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 found a solution to this problem.
After having installed SRA toolkit:
prefetch -v SRR9....1
this helps download the SRA file.
Then to extract it as FASTAQ
fasterq-dump SRR9....1

That just postpones the problem. Now you need to align these files, tell me a single popular Widows aligner.

@ATpoint clustal? Please try and make your comments helpful instead of just snide. Of course you cannot seriously do bioinformatics if you just use Windows, but people who only dabble might be able to get away with it, and the specific example of alignment tools is a bad one.


text:
Fetching uniprot seq using bash command
I am trying to get a protein sequence using this command:
curl https://www.uniprot.org/uniprot/A2Z669.fasta

But I don't get any output:
$ curl https://www.uniprot.org/uniprot/A2Z669.fasta
$

And no file is created.

The URL https://www.uniprot.org/uniprot/A2Z669.fasta redirects to https://rest.uniprot.org/uniprotkb/A2Z669.fasta. You can see this if you open it in your browser: clicking on the first link will take you to the address of the second.
Now, curl doesn't follow redirects by default, you need the -L option for that (see man curl):

-L, --location
(HTTP)  If  the  server  reports  that  the requested page has moved to a different location (indicated with a Location:
header and a 3XX response code), this option makes curl redo the request on the new place. If  used  together  with  -i,
--include or -I, --head, headers from all requested pages are shown.

So you need -L for curl to work here:
$ curl -L https://www.uniprot.org/uniprot/A2Z669.fasta
>sp|A2Z669|CSPLT_ORYSI CASP-like protein 5A2 OS=Oryza sativa subsp. indica OX=39946 GN=OsI_33147 PE=3 SV=1
MRASRPVVHPVEAPPPAALAVAAAAVAVEAGVGAGGGAAAHGGENAQPRGVRMKDPPGAP
GTPGGLGLRLVQAFFAAAALAVMASTDDFPSVSAFCYLVAAAILQCLWSLSLAVVDIYAL
LVKRSLRNPQAVCIFTIGDGITGTLTLGAACASAGITVLIGNDLNICANNHCASFETATA
MAFISWFALAPSCVLNFWSMASR

Or, alternatively, you can use wget which does follow redirects by default:
$ wget "https://www.uniprot.org/uniprot/A2Z669.fasta"
--2023-12-10 16:21:46--  https://www.uniprot.org/uniprot/A2Z669.fasta
Loaded CA certificate '/etc/ssl/certs/ca-certificates.crt'
Resolving www.uniprot.org (www.uniprot.org)... 193.62.193.81
Connecting to www.uniprot.org (www.uniprot.org)|193.62.193.81|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://rest.uniprot.org/uniprot/A2Z669.fasta [following]
--2023-12-10 16:21:46--  https://rest.uniprot.org/uniprot/A2Z669.fasta
Resolving rest.uniprot.org (rest.uniprot.org)... 193.62.193.81
Connecting to rest.uniprot.org (rest.uniprot.org)|193.62.193.81|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://rest.uniprot.org/uniprotkb/A2Z669.fasta [following]
--2023-12-10 16:21:46--  https://rest.uniprot.org/uniprotkb/A2Z669.fasta
Reusing existing connection to rest.uniprot.org:443.
HTTP request sent, awaiting response... 200 OK
Length: 314 [text/plain]
Saving to: ‘A2Z669.fasta.1’

A2Z669.fasta             100%[=======================================>]     314  --.-KB/s    in 0s

2023-12-10 16:21:47 (8.57 MB/s) - ‘A2Z669.fasta.1’ saved [314/314]

That creates a A2Z669.fasta file in the current directory which will contain the sequence.


text:
Fst calculation from VCF files
I have four vcf files, SNPs_s1.vcf, SNPs_s2.vcf, SNPs_s3.vcf, and SNPs_s4.vcf, which contain information about SNPs. These vcf files were obtained by using the following methods: the initial input files were short-paired reads

I did mapping with minimap2
./minimap2 -ax sr ref.fa read1.fq.gz read2.fq.gz > aln.sam

converted to bam file
samtools view -Sb aln.sam > aln.bam

Sorted the bam file
samtools sort aln.bam > aln_sorted.bam

created an index file
samtools index aln_sorted.bam

created an index file of reference
samtools faidx ref.fa

extracted the raw variants using freebayes
freebayes -f ref.fa aln_sorted.bam > out_var_s1.vcf

extracted SNPs using vcf tools
vcftools --vcf out_var_s1.vcf --remove-indels --recode \
   --recode-INFO-all --out SNPs_s1.vcf

I would like to calculate the fixation index (Fst) between SNPs_s1.vcf and SNPs_s2.vcf, SNPs_s1.vcf and SNPs_s3.vcf etc. I attempted to calculate Fst using bcftools, but I have to provide the id list. I tried the code below for that.
bcftools query -l snps_s1.vcf

The output was "unknown".
How do I calculate Fst using the above vcf files? Your suggestions would be appreciated.

Can you please explain a bit more about Fst, and why you need to calculate it? I'm used to Fst being interpreted as a population statistic, and it doesn't [currently] make sense to me to calculate it for only a few individuals.

@gringer I would like to use Fst to find the genetic diversity between populations. Which method is better for few individuals?

Who is asking you to calculate Fst, and why? If you're doing population genetics, you need populations. A few individuals do not make a population, unless it's kākāpō or kakī, or some other similarly endangered species.

The first issue is easy, the name of your sample in your snps_s1.vcf is set to unknown. You can fix that bit with sed -i '/^#CHROM/s/unknown/sample1/' snps_s1.vcf. But gringer's points are the main issue. What are you really trying to do here with only 4 individuals?


text:
Using Diffusion Map on normalized - or PCA reduced data?
This question was also asked on GitHub
In the R package destiny, a Diffusion Map is often performed on the first 50 PCs of the data.

Is there a theoretical reason for performing a Diffusion Map on PCs rather than on the data itself ( other than speed)?

Is it reasonable to do the diffusion map on the normalized data itself rather than on Principal Components of the data, if that differentiates the states better ?

Techinically, I am using the  R Bioconductor package destiny. The help page on destiny::DiffusionMap tells "Using e.g. 50 PCs results in more regular looking diffusion maps", but it's not clear compared to what will they be more regularly looking.

The empirical statistical point of view
The technique is  uniform manifold approximation and projection (UMAP) and is comparable to t-SNE. UMAP is considered a slight improvement, personally I think they have similar behaviour. They are dimensionality reduction methods.
The reason UMAP or t-SNE is often not applied directly is it is weaker than PCA are capturing the multi-dimensional variance of the data. Both methods can be applied directly to the data, but will not necessarily deliver stronger clusters than PCA (which is the goal of the technique).
However, the weakness of PCA is that the variance can be spread across 50 principle components, making it hard to recover the representative clusters in certain situations. These situations are commonly found in RNA expression analytics, notably where the first 3 or 4 PCs only comprise ~30% of the variance. PCA should be trapping >80% of the variance within the first 3 or 4 PCs. In expression data (RNAseq) the next 46 PCs might comprise a total of ~30% of the variance and omitting that is mistake, but manually comparing 50 PCs in bivariate plots is not practical. So at best a good unsupervised learning algorithm (UMAP/t-SNE) following PCA in RNA expression via RNAseq might trap 60% of the variance. It's not ideal but it's okay. Thus you can't ignore the 30% of the variance in the bottom 46 PCs and UMAP/t-SNE is recovering that. If 80% of the variance was in the top 4 PCs, then the remaining 20% of the variance in the lower 46 PCs can be ignored and t-SNE or UMAP is not needed.
Resolution The way you can capture the variance from all 50 PCs is to perform a clustering method over the top of PCA, this being t-SNE or UMAP.
Summary,

t-SNE or UMAP without PCA might trap at best 20% of the variance.
The first 4 PCs in a PCA might trap 30% of the variance (it will be a better result than t-SNE or UMAP).
Combining PCA and then UMAP/t-SNE will capture 60% of the variance.

60% of the variance, its just about acceptable. Thus thats why this approach  is done from an empirical point of view.
Practical point of view
PCA and then UMAP (or t-SNE) will give tighter, neater clusters (cleaner, less messy groups), but thats not the reason why you are using it.

The publications in the link,
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10312943/
and the PLoS Biology paper
https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001532
... both use UMAP across PCA. Its a major approach to cluster resolution.
Diffusion maps are an attempt to extrapolate a variance within Eigenvalues via a random walk, i.e. to defuse within Euclidean space. It is not a primary calculation in my personal opinion because you are sort of introducing more variance into a variance dimensionality reduction calculation. Once you have very tight clusters using PCA/UMAP or PCA/tSNE then you could check to see whether Diffusion map clusters are similar.
This is my perspective only, not to replace PCA with Diffusion map nor to consider a "down-stream" calculation. Its separate calculation used to verify if PCA-UMAP is okay.
If you doing graph theory a Diffusion map could be really important as a null hypothesis and a really important way to obtain a different understanding of the graph (i.e. this is a technical term for a network). Using a diffusion map without any directionality (the papers you cited)*, I struggle to see this as a primary analysis. In the PLOS Biology paper is part of the appendix.
*, graph says A->B->C
The publications treat A,B, C as unconnected. If you can hypotheses how A, B and C are connected then diffusion map could be cool. Thats my opinion from the perspective of graph theory.
Summary I would definitely choose PCA-UMAP over a diffusion map every time. This is what has happened in the two publications cited. The diffusion map is secondary analysis, so by all means do it ... but present PCA-UMAP first, and base the results around that.

From the comments there's some sort of time-series. So if you are performing a time-series - it's not datamining - this would be the justification of using a Diffusion map over PCA. This therefore answers points 1 and 2. However, it would imply you are pooling your data.

The specific technique is a Diffusion Map after PCA (and not UMAP/t-SNE). Does the same reasoning hold for the Diffusion Map dimensionality reduction as well?

Thanks for clarification, I've answered above

I'll explain my specific case in more detail : The reason I personally use a Diffusion Map is not only to see the clusters, but to draw a trajectory between the clusters. In my case, not only the Diffusion Map separates better between the (two) main clusters, but the genes change in a meaningful way, across the trajectory : from state A, to intermediate state, to state B.  However, this is not stable in the sense that it appears only when using the normalized data as a basis for Diffusion Map, not when using the 50 PCs.Hence, the two OP questions focus on Diffusion Map(on normalized vs on PCs)

@Sam maybe this is a separate question? The studies you presented suggested a straight cluster analysis. If you are doing e.g. a time series that will require a lot more thought.

I have intended the studies just to illustrate the point that taking the first 50 PCs is often done (so that I do not get called out for making unfounded statements about what is often done). They were not intended to be the focus of the question, but I see it was not clear.
How would you phrase a separate question? In my mind, the question is well-phrased (but it might be I don't get it).

Are you doing a time-series?

My data has been collected at one time-point, but computationally I perform trajectory analysis using Slingshot, based on that, divide the trajectory into 3 timepoints (start, mid, end), and perform differential expression test between the timepoints using tradeSeq

can you help further?

Thanks @Sam answered above


text:
What is the troubleshoot for this error: conversion of .SRA to FASTA file on command prompt?
I am getting this error message after using the following code:
C:\sratoolkit.3.0.7-win64\sratoolkit.3.0.7-win64\bin>fastq-dump --fasta SRR1658345

Error:
2023-12-11T06:08:04 fastq-dump.3.0.7 err: timeout exhausted while waiting condition within process system module - failed SRR1658345

=============================================================

An error occurred during processing.
A report was generated into the file 'C:\Users\Hp/ncbi_error_report.txt'.
If the problem persists, you may consider sending the file
to 'sra-tools@ncbi.nlm.nih.gov' for assistance.

=============================================================

2023-12-11T06:09:46 fastq-dump.3.0.7 err: connection failed while opening file within cryptographic module - Cannot KStblHttpFileTimedReadChunked: to=0

What should be done?
Link to ncbi error report
Edit:
I am on Windows, using Command Prompt.

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.

Please [edit] your question and explain what you are doing. First of all , tel us your operating system. Are you running Windows or some reason? If yes, what exactly are you running? Is that WSL? Or are you running it from a Windows shell? Also, please upload that error report to something like https://pastebin.com instead of making it a pdf. We need to be able to copy text from it and that site doesn't allow copying.

Isn't the syntax    --fasta <line-width> where <line-width> can be                                      "default" or "0" for no wrapping"

@Sam I don't think so. I just tried running fastq-dump --fasta SRR1658345 on my Arch Linux machine and it brought me a fasta file just fine.


text:
Uploading manifest from CCDI to CGC project
Thank you in advance for your help!
I am trying to download data from the CCDI (Children Cancer Data Initiative) and I have a manifest csv file that I downloaded from the CCDI. I'm trying to upload this data to my CGC project, but there is not option to import data from manifest file from the CCDI specifically, and when I do try to upload it, only the manifest csv is added and not the files. What is the process for uploading manifest that are not from one of these databases?

Thank you so much!

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.

Please [edit] your question and give us (a lot) more context. What is this screenshot from? Where could we go to find it? Why are you using "import", is there no "upload" option? Do you have a manifest file or a csv file? Can we see an example? What is CGC? Do you mean https://cgc.umn.edu/? Or maybe https://cancergenomics.org/? What would uploading these files actually do? What kind of data do you have?


text:
How to remove third codon positions from a charset in iqtree?
I need to build a phylogenetic tree using IQ-TREE, starting from a sequence alignment in CODON format of several invertebrate mitochondrial genes. These are my charsets:
nexus
BEGIN SETS;

charset atp6.fasta = 1-651;
charset atp8.fasta = 652-729;
charset cob.fasta = 730-1749;
charset cox1.fasta = 1750-3237;
charset cox2.fasta = 3238-3858;
charset cox3.fasta = 3859-4632;
charset nad1.fasta = 4633-5484;
charset nad2.fasta = 5485-6156;
charset nad3.fasta = 6157-6459;
charset nad4.fasta = 6460-7683;
charset nad4L.fasta = 7684-7899;
charset nad5.fasta = 7900-9294;
charset nad6.fasta = 9295-9663;
END;

I would need to remove the third bases, how should I write them? These are the commands with which I would launch IQ-TREE:
iqtree -s "fasta_file" -bb 1000 -alrt 1000 -m MFP+MERGE -spp "charset.nexus" -st CODON5

Please, any suggestions would be very helpful.

Yes, it's easy.
#nexus
BEGIN SETS;

charset atp6.fasta = 1-651\3, 2-651\3;
charset atp8.fasta = 652-729\3, 653-729\3;
charset cob.fasta = 730-1749\3, 731-1749\3;
# ... so on for every gene
END;

So for atp6 you will get 1,2,4,5,7,8,10,11 .... up to 651, and so on for each gene. Thus the residues at position 3, 6, 9, 12 are excluded.

Thank you very much, but by using the above charsets I get  "ERROR: Range 1-651\3, 2-651\3 is not multiple of 3 (necessary for codon data)". It's the same error that I get using other solutions suggested on the web; I don't understand why, using regular charsets without removing the third bases everything runs smoothly.


text:
Why are Minimap2 alignments different with CIGAR generation flag?
I am using Minimap2 (v2.26-r1175) in Linux to generate a sequence alignment between the Streptomyces coelicolor A3(2) chromosome (ref.fa) and the Mycobacterium tuberculosis chromosome (query.fa). My desired output is a PAF (Pairwise mApping Format) file.
The general way to align reference and query sequences with Minimap2 is the following:
minimap2 ref.fa query.fa > approx_mapping.paf
# OR equivalently...
minimap2 -k15 -w10 ref.fa query.fa > approx_mapping.paf

You can get Minimap2 to generate custom CIGAR tags (cg:) in the PAF by adding the -c flag:
minimap2 -c ref.fa query.fa > approx_mapping_cigar.paf
# OR equivalently...
minimap2 -c -k15 -w10 ref.fa query.fa > approx_mapping_cigar.paf

Moving to an R environment [with the intention of downstream data visualization], you can quickly assess the number of observations in each object (a.k.a. rows in each dataframe) with dim():
file1 <- "approx_mapping.paf"
paf_basic <- read.table(file1, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file1, sep = "\t")))))
dim(paf_basic)
[1] 205  18

file2 <- "approx_mapping_cigar.paf"
paf_cigar <- read.table(file2, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file2, sep = "\t")))))
dim(paf_cigar)
[1] 200  24

Note the unqiue way we have to read the PAF files into R using read.table() to ensure no columns are cut off and added to new rows. This is because read.table() is designed to examine only the first five lines of an input file to determine the number of columns in the dataframe. Any columns on lines beyond the fifth line that exceed the initially determined count are appended to a new row. In the case of PAF files generated by Minimap2 with CIGAR tags, the number of columns can vary per line, which makes this necessary.
With the default settings, I get 205 observations and with the addition of the CIGAR flag, I get 200 observations.
Why do I get fewer computed alignments when the CIGAR tags are added? I would expect the alignments between the reference and query to be identical between these runs.

Comparing the PAF files (directly generated by minimap2) after sorting (sort -k3,3 -k4,4) using vimdiff reveals that the alignments computed between the default and CIGAR setting runs are largely similar, but not identical:

Note that the image above only shows the first 44 rows of both files.
Due to the small differences seen in alignment coordinates between these two dataframes (see fields $3 and $4 in the vimdiff screenshot), it is NOT trivial to determine which five observations (rows) are new alignments computed when using the CIGAR flag -c.

Plotting these dataframes confirms that the paf_basic and paf_cigar objects do, in fact, contain some different sequence alignments from each other:
library("ggplot2")
ggplot(data = paf) +
  geom_segment(aes(x = V8, xend = V9, y = V3, yend = V4), size = 2) +
  labs(title = "paf_name", x = "Reference (Mb)", y = "Query (Mb)") +
  coord_fixed(ratio = 1) +
  theme_bw()

So, it appears as though minimap2 -c produces a slightly different alignment than minimap2, which is surprising to me and begs the question: Which alignment is better or preferable?

You've shown that you have 5 extra rows when you read it into R. However, you haven't actually shown what the extra rows are. If we could see them, we might be able to guess what's going on. I can't reproduce the issue on my own. Can you rule out the trivial explanation that R's read.table is parsing the PAF incorrectly? That seems most likely given your observation about tag-only rows.

Great point! Unfortunately, there is no straightforward way to determine what the 5 extra rows are because the two data frames have no matching observations between them (checked by dplyr::semi_join(paf_basic, paf_cigar, by=c("V3","V4","V8","V9"))). Based on the graphs included in the question (a recent edit), it seems that many of the alignment coordinates are close, but not identical. Unfortunately, the dput() outputs for these data frames are too large for SO question character limits. However, I have verified that the *.paf files match the paf R objects produced by read.table().

have you looked at the two PAF files side by side in e.g. vimdiff? It should not be too hard to find the differences leading to 2 extra rows. But IDK exactly what is going on. If you've demonstrated that adding the -c option changes the alignments themselves,  that certainly sounds wrong. I'd suggest posting an issue on minimap2 github, Heng is pretty responsive. (If you do so, be sure to link to the issue here!)

Thanks for the idea Maximilian -- I hadn't heard of vimdiff before. I gave it a go and found it immensely useful. I've updated the question to include a snapshot of the two files sorted and side-by-side in the program which reveals just how similar (but not identical) the data frames are. I'll look into posting the issue on the minimap2 github. Thanks for taking time out of your day to ponder this issue and help troubleshoot with me. Cheers!

Note also that the scores are very different with and without -c. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea

Awesome analysis Colin! Interesting to see how much of an effect -c has on identity scores of different length alignments. Thanks for adding it to the conversation here.

The program's author pointed me to the Minimap2 FAQ:

1. Alignment different with option -a or -c?
Without -a, -c or --cs, minimap2 only finds approximate
mapping locations without detailed base alignment. In particular, the
start and end positions of the alignment are impricise. With one of
those options, minimap2 will perform base alignment, which is
generally more accurate but is much slower.

So, to specifically address the main questions asked here:

Why are there fewer and non-identical computed alignments when the CIGAR tags are added?

minimap2 alignment accuracy differs when running minimap2 -k15 -w10 vs. when running minimap2 -c -k15 -w10. It is actually not unexpected that the alignment coordinates in each PAF file are only similar, not identical, and an increase or decrease in the number of alignments between the outputs is simply an artefact of this accuracy discrepancy.

Which alignment (default vs. with CIGAR tags) is better or preferable?

Although both alignments report very similar sequence relationships, adding the -c flag for CIGAR tag generation is preferable to get more accurate base alignments.

Thanks for posting the answer and accepting it. Kudos to you @Gawain

Note also that the "score" can vary dramatically depending on whether -c is used or not. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea it is also alluded to briefly in the PAF.md where score is a required field but "If alignment is not available, column 10 and 11 are still required but may be highly inaccurate."


text:
What method uses "diverge probability" in RNA-Seq differential expression analysis?
I am replicating differential expression analysis from a paper (here) on RNA-Seq data. I extracted the regulated genes using DESeq2::results with a threshold of 1 for log2 fold change and 0.05 for the p-value. The results are very different from that of the paper. The paper does not mention much information about their statistical methods except for this paragraph:

Finally, the libraries were sequenced using the BGISEQ-500 platform. Differentially expressed genes were identified based on fold change ≥ 2 and diverge probability ≥ 0.8. The functional categories of the differently expressed genes and the pathway were obtained using MapMan as the classification source.

Now, I have no idea what they mean by "diverge probability" but it looks like another method to get the regulated genes. My questions are:
Question 1: What could be the method they are using in their analysis?
Question 2: Are the results from the "diverge probability" method different from the ones we get using DESeq2::results?
Question 3: If the answer to Question 2 is yes, how do we get from one to the other?
I appreciate any insights or experiences you can share regarding this approach. Thank you!

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'm pretty informed about stats but haven't a clue what "diverge probabilities" are. Are you referring to Kullback–Leibler divergence by any chance? Whilst I understand the question, at present its unanswerable but the metric doesn't appear to correspond to a known critical value approach.

@M__ I tried to add more info and edit my questions to be more clear. Hope that helps, and let me know if you want me to add more info.


text:
PCR primers for detecting hetero vs homozygous and WT mice in our colony
My PI has asked me for advice on making a PCR screen to identify mice from a particular Cre line.  We have multiple Cre lines in the colony so it can't just be for mice with a copy of Cre in their genome the primer. The insertion site is known for the mutation,  is there a way to use blast or some other search tool to pull the sequences we should base our primers on?
We also want to identify which mice are hetero vs homozygous.  Are there particular kinds of PCR we should use or will a well designed primer enable this?

Details
The strain:  https://www.jax.org/strain/004146
I was thinking of using this tool:  https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome
The strain reference from Jax, not sure the recommended cre primers are going to be sufficient for what we want.  So I'm trying to find the FASTA file for this strain in the region of interest, then I can throw that into the BLAST tool for primer design I listed above.  Not sure where to search for the FASTA for this strain though. Strain information: https://www.informatics.jax.org/reference/J:66884
Jax recommends this protocol so I suggested it:  https://www.jax.org/Protocol?stockNumber=004146&protocolID=20627

Jax recommends this protocol so I suggested it:  https://www.jax.org/Protocol?stockNumber=004146&protocolID=20627

Primer-blast is okay, but I presume these are lab mice so the variation is low, thus  you can just use Primer3. SNPs are detected via TaqMan here, i.e. a separate probe that sits over the mutation site, it will detect heterozygosity via qPCR signal (half strength than homozygous).


text:
Making a matrix when samples that have missing data
I am trying to make a dataset for methylation analysis.
I called methylation using Methyldackel and there are various positions without methylation calls. This creates a problem because when I try to use for example, DMCHMM or MethylSig, to estimate differentially methylated regions, the initial step prompts an error:
Building BSData object.
Warning message:
In .merge_two_Seqinfo_objects(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': NW_022883351.1, NW_022883409.1
  - in 'y': NW_022883010.1
  Make sure to always combine/compare objects based on the same reference
  genome (use suppressWarnings() to suppress this warning).

I used the same reference for all my samples, but when it tries to combine the files, it complains because some samples don't have data (missing) in some positions.
Is there a way to use R or Linux to fix this, and check those lines positions with missing data and add those (lines) positions in the files?
The files look like this:
chrBase chr     base    strand  coverage        freqC   freqT
NW_022882922.1.4856     NW_022882922.1  4856    F       2         0.00  100.00
NW_022882922.1.4866     NW_022882922.1  4866    F       2         0.00  100.00
NW_022882922.1.4889     NW_022882922.1  4889    F       4         0.00  100.00
NW_022882922.1.4892     NW_022882922.1  4892    F       4         0.00  100.00

so that all the files have "information" in all the positions?
Thank you
Juan Pablo


text:
align to the whole hg38 genome, then split bam and collect metrics on each bam issue
I am doing some coverage analysis on deep sequencing human genome. I first align to the whole hg38 genome, then split bam to each chrom,  and collect metrics on bam for each chom.
I first split using samtools. and then tried bamtools.
samtools: samtools view -bo ${pair_id}.\\\$chr.bam $bam \\\$chr
bamtools:bamtools split -in $bam -reference

When I collect metrics I am using the reference fasta from each chrom. The picard metrics complains the index dictionary size is different:  Sequence dictionaries are not the same size (25, 1). The picard method I used is: MarkDuplicates, CollectAlignmentSummaryMetrics, CollectGcBiasMetrics, CollectInsertSizeMetrics, CollectWgsMetrics.
My data is really huge, 70G for a sample, that is why I would like to split when doing the picard metrics. Seems picard metrics collection use a lot of memory. So for the 70G human fastq(30X), how big the memory is needed for running picard wgs metrics?
How can I deal with this? Can I modify the split bam header? How to do it?
Or not split it into chromosmes, rather, using an interval list when collecting the metrics? will that use less memory?

Please [edit] your question and show us the exact commands you are using and the exact error message. My guess is that you're giving the index for the entire genome, so the error makes sense, but we can't know without the commands. Can you explain why you are splitting the bam though? Why not get the stats for the entire bam, and then split those by chromosome?

Can you please explain what metrics you are interested in? If it's something equivalent to "number of reads mapped per chromosome", there's probably a better way to get those statistics.

Thanks for the edit, but we still don't know what Picard commands you are giving. What index files are you using? Is it an index for the entire genome and you're feeding only a per-chr fasta file? Are you referring to the bam index? Have you indexed the individual bam files? Also, 70G isn't huge, that's a pretty normal file size for this kind of thing, so you might want to consider more powerful hardware if that is an option.

Thank you! So first I aligned with the whole genome, and then tried to split the bam into different chromsome, then to do the picard using the per chrom ref.  Thank you for your comment on the size(70G not huge), I have decided to just to use large memory rather than split them. I am OK if this post is to removed from this site, since it might not be helpful for anyone

This can be done, but you have to split the BAMs while keeping the full header (with the dictionary) in place for all the shards. This can be done with gatk PrintReads. One thing to remember is that you all need to have one shard with all the "unmapped" reads. Otherwise you will not have the complete information. You can get the unmapped reads by putting unmapped into the -L argument.
That said, MarkDuplicates cannot be run on shards since it needs both reads from every pair and some of the read pairs might be split. So you'll need to mark-duplicates on the entire 70GB bam (sorry).
You could split the BAM into shards but only the "proper pairs" ( filter after GATK with samtools view, but only on the sam flag so you don't mess with the dictionary) which will make sure that your pairs are on the same contig. Then put all the not proper pairs into their own shard.
You can also do the proper/improper filtering with GATK, but it involves read filters which are a bit more complex than samtools flag filters.


text:
How can I use QIAGEN Ingenuity Pathway Analysis (IPA) on proteomics data associated with an outcome?
I am trying to do pathway analysis on a set of proteins associated with several outcomes (one is continuous, the others four are binary). I have run independent associations between the proteomics data and the outcomes in question and have ascertained betas, hazard ratios, p values, and associated multiple-testing correction p-values. Out of the set of 3000 proteins for which we have proteomics data, we have found significant FDR corrected associations with our outcomes in ~1800. I am trying to do pathway analysis with these proteins.
I am trying several pathway analysis approaches, including fgsea, online GSEA through GenePattern, and Ingenuity Pathway Analysis. For the GSEA approaches, I have been able to generate the ranked list and provide the ranked list to the software. However, it is my understanding that IPA is a form of overrepresentation analysis and therefore I should not provide a ranked list. However, when uploading my dataset to the IPA application, it asks for expression change values (my understanding is that all these approaches were originally designed for microarray-based  differential expression data).
Since we are looking at the association between proteomics and several outcomes, expression change values are not relevant to these data. However, the directionality of the effect is relevant when considering pathways, and since we have the betas for the association it seems to me that I should include it. I looked through the IPA documentation and all of it seems specific to differential expression analysis. Can anyone point me in the correct direction to perform IPA on the association between my outcomes and proteomics data? Any help would be appreciated.

What have you already asked Qiagen about this? IPA is a very expensive product; Qiagen should at least provide good support around the bioinformatics side of things because bioinformatics is such a core part of their product.


text:
tBLASTn vs NCBI TSA
Is there a way to blast my proteins against all the TSA (Transcriptome Shotgun Assembly Sequence Database) belonging to a specific taxonomy using command line?
Thanks,
Marco

This seems like an XY problem. Please add more detail to your question post, including information about the problem that you are trying to solve by doing this. Have you considered other approaches to discover the taxonomy of proteins?

Just go to Blast here

Click where is says Database go to Transcriptome Shotgun Assembly Sequence Database (TSA).

Click where is says "Limited by" click organism and thereon specify the group of interest.

Thats it.

The OP wants to do this "using command line" though. And I haven't managed to find a good way to do that, personally, which is why I came up with that really complex workaround for a similar question.


text:
Prediction of DNA binding motif from sequence/structure of TF/regulator protein?
Do you know if it is actually possible to predict with good accuracy a DNA binding motif from only the sequence/structure of the DNA-binding protein?

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.

This question looks like it's too general to be useful for asking or answering; please frame your question around a specific problem that you're having, especially including detail about the nature of the DNA binding motif that you're interested in.

Find papers....

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


text:
Volcano Plot: How to show unique gene names in the plot?
I'm wondering how can I show unique gene names in the volcano plot.
My data looks like,
head(DMCf1_huMetGNc)
     chr    start    end.x       pvalue       qvalue meth.diff        gene
152 chr1 24909805 24909805 1.915400e-03 0.0168087962 -52.00000 RUNX3;RUNX3
168 chr1 24931440 24931440 8.645625e-03 0.0489503201  28.98551 RUNX3;RUNX3
179 chr1 26862096 26862096 1.308307e-04 0.0023224596 -42.48366         SFN
187 chr1 26914178 26914178 4.508416e-06 0.0001786963 -57.97872       NR0B2
188 chr1 26914203 26914203 5.337425e-05 0.0011910452 -47.30424       NR0B2
194 chr1 35573100 35573100 1.433689e-03 0.0135780139 -42.04545      TFAP2E

The original data file is quite large, here under the gene column some of the rows have multiple same gene names separated by ;. How can I show only unique gene names.
No pattern of gene names exists when multiple names are present. In certain cases, the same gene is repeated, while in other cases, distinct genes are present in the same cell. The majority of the time, a single gene exists in a cell.
library(ggplot2)
library(ggrepel)
    v1<-ggplot(data=DMCf1_huMetGNc, aes(x=meth.diff, y=-log10(pvalue) )) + geom_point()
    v1+geom_text_repel(aes(label= gene))

Here, you may see how the FGFR2 gene is repetitive in the plot that I want to remove and the same name from a single row should appear only once in the plot.

What do you want done with lines that include multiple gene names?

I want to show the unique names only, not the same repetitive  names.

So if there's A;B;A, you wish to show A;B and if there's A;A, you want just A?

This explanation is not descriptive enough. Do you want lines to be duplicated, one for each gene? Do you want to exclude lines with multiple genes? Or something else?

Yes,  RamRS is true. I want to show the names as @RamRS is suggested.

Based on your comments, you wish to obtain a semi-colon separated list of unique values in the gene column. You can create a new column containing these using:
xyz <- data.frame(chr="chr1", gene = c("RUNX3;RUNX3", "RUNX3;RUNX3", "SFN", "NR0B2", "NR0B2", "TFAP2E", "A;B;A"))

xyz
   chr        gene
1 chr1 RUNX3;RUNX3
2 chr1 RUNX3;RUNX3
3 chr1         SFN
4 chr1       NR0B2
5 chr1       NR0B2
6 chr1      TFAP2E
7 chr1       A;B;A

xyz$gene_uniq <- sapply(xyz$gene, function(x) paste(unique(unlist(strsplit(x, ";"))), collapse=";"))

xyz
   chr        gene gene_uniq
1 chr1 RUNX3;RUNX3     RUNX3
2 chr1 RUNX3;RUNX3     RUNX3
3 chr1         SFN       SFN
4 chr1       NR0B2     NR0B2
5 chr1       NR0B2     NR0B2
6 chr1      TFAP2E    TFAP2E
7 chr1       A;B;A       A;B

Or you can use dplyr. I'm adding a column called entry_num below since my example has non-unique combinations of chr and gene that we need to retain.
library(dplyr)
library(tidyr)

xyz <- data.frame(chr="chr1", gene = c("RUNX3;RUNX3", "RUNX3;RUNX3", "SFN", "NR0B2", "NR0B2", "TFAP2E", "A;B;A"))
xyz %>% mutate(entry_num = 1:nrow(xyz)) %>% separate_rows(gene, sep = ";") %>% unique() %>% group_by(chr, entry_num) %>% summarize(gene = paste(gene, collapse=";")) %>% select(-entry_num) %>% as.data.frame()

   chr   gene
1 chr1  RUNX3
2 chr1  RUNX3
3 chr1    SFN
4 chr1  NR0B2
5 chr1  NR0B2
6 chr1 TFAP2E
7 chr1    A;B


text:
pacbio demultiplexing issue
I was wondering if there is a way to provide pacbio lima demultiplexing only a barcode, and let  lima demux sample even if only one end has the barcode, not both?

Can you please explain more about the context of your problem? What are you trying to do with the demultiplexing; what is the ultimate aim of the analysis?

Some examples of what your data look like might help

PacBio lima need both 5'end and 3' end barcodes. You can have a unique barcode in one end and multiple on the other end; but it need two barcodes.


text:
Longstitch error make: command: Command not found *** No rule to make target
I installed Longstitch and ran the test script with no issues. The output files matched the expected output files. But when I am now trying to run Longstitch on my own data I am getting this error.
(medaka) operations@LabPC1:~/medaka/longstitch-1.0.5$ ./longstitch run draft=~/medaka/~medaka/medaka_results/consensus.fa reads=~/medaka/combined.fastq G=4.5e8
make: command: Command not found
make: command: Command not found
make: *** No rule to make target'/home/operations/medaka/~medaka/medaka_results/consensus.fa.fa', needed by'/home/operations/medaka/~medaka/medaka_results/consensus.fa.cut250.tigmint.fa'.  Stop.

I changed the file name from .fasta to .fa for my draft assembly because it requires a .fa file, but now I am getting ".fa.fa" What is wrong here?

Please provide a link and/or additional information about what Longstitch is. It would also be useful if you could provide a few lines of input and output (e.g. from both the example data and your own data).

Could you download it from conda and try again? Its here, https://docs.conda.io/en/main/miniconda.html
Then,
conda create --name myenv
conda activate myenv
conda install -c bioconda longstitch
longstitch ...

Conda avoids all the build stuff.

Notes
If you find that "conda" is slow and doesn't ever load.
There is an alternative to conda called mamba, this is fast and designed to address this exact issue. I have never needed to use it so the following is a rough guide only.
Download and install via
conda create --name myenv
conda activate myenv
conda install -c conda-forge mamba

The documentation is here, https://mamba.readthedocs.io/en/latest/user_guide/mamba.html
I've never used it but the idea is ...
mamba create -n mambaenv <list of packages>
mamba activate mambaenv

The list of packages could be left blank for now, then,
mamba install "Python>=3.10.0" numpy pandas pathlib
mamba install longstitch

Can also do
mamba install -n mambaenv longstitch

NOTE longstitch exists strictly in -c bioconda. You do not seem to search and install via channels in mamba
You search for a package and its dependencies via
mamba repoquery search longstitch
mamba repoquery depends -a longstitch

It is important to place a virtual environment between the system Python and package Python (don't want a system corruption).

Finally ...
For fa.fa ... hmmmm, try ls *.fa and check for .fa.fa files.

The main issue I can see is the filenames, and should be easily solvable.
Simply put, LongStitch is very particular about filenames. It expects the input draft to have the extension .fa, and the input reads to have any from a range of extensions (your .fastq is perfectly fine, though).
The main issue is that LongStitch expects the extensions to be missing from the command. So if your draft is draft is called consensus.fa and your reads are  combined.fastq, then the command would include something like draft=consensus reads=combined, not draft=consensus.fa reads=combined.fastq.
(unsure if it accepts input from different directories, or whether you'll need to move your input to the same directory as well)
More generally, when having issues with a particular tool which is available via github, I would encourage you to have a look at the 'Issues' or 'Help' section. What I normally do is to delete the 'is:issue is:open' terms in the search bar, which allows me to search for my bug among older, potentially solved issues (e.g. https://github.com/bcgsc/LongStitch/issues/54), and also generally to check that the software is still supported.
Good luck!
(*that said, M_'s suggestion of using mamba is generally a good one, and has helped me with another LongStitch bug)


text:
Snakemake and github actions
I am working on a snake-make workflow that follows the guidelines of formatting and linting so that it can be published on the official Snake-make repository.
Whenever I push my code on GitHub, the GitHub actions are run and the pipeline is tested as indicated here.
The problem is related to the execution of some time-consuming rules, such as the genome indexing generation for bowtie. These rules are run every time and take plenty of time.
I was wondering whether there was a way to bypass the execution of these rules and test the pipeline anyway. Do you have any suggestions about this?
A problem is that I have a rule in my snakemake workflow that downloads the human genome and that is needed to index and then for the alignment. I need a workaround to avoid this behaviour, since has you said it is not feasible to test on the whole genome.
Thanks

Please add more context to your question. An example snakemake workflow added into your question with code formatting would be helpful (i.e. not just linking to an external site).

The answer for any testing purpose is generally to put together a small and minimal dataset so the test completes in seconds or minutes while still being representative. It makes (to me) limited sense to have a Actions/CI running using for example the entire human genome as input. The point of testing is to check that steps are executed properly, arguments are passed correctly, errors are detected. For this use either a small genome (such as from a bacterium or yeast), or make an artificial one consisting of very few basepairs.
What I often do (for example for a peak calling workflow) is to take a full dataset, then subset to a region where I see a) many peaks and b) no peaks, then extract this genomic range (a few thousand bp) to build an artificial genome, plus extract from the bam/fastq the reads contributing to this region. That way I have basically a read-world dataset but very limited in size. For CI testing that does very well and CI finished in seconds, allowing to test even different configurations of this workflow on a standard / unpaid Actions runner.
I personally prefer to not automate downloads but rather to have the user provide an actual file, or provide the file pre-downloaded in the workflow bundle. Connections might be bad or broken that day or the server/HPC is not connected to internet, and the entire workflows fails. Maybe write the relevant section in a way that it either accepts a URL for download OR a local file, and then write a testing profile that then reads the file to ensure it always is accessable to the Actions runner.


text:
otu_table from phyloseq object does not have correct abundance: all asv are "0" for specific group of samples
This question was also asked on stack overflow
I have QIIME2 outputs (feature table, tree, taxonomy, and metadata) that I've imported into R using qiime2R (v. 0.99.6) ‘qza_to_phyloseq’ to create a phyloseq object. My QIIME feature table shows reasonable values for ASV abundance across all samples (n=16; feature table available upon request).
Problem: The otu_table the ASV values for an entire study (NCBI SRA accession PRJNA217520, named AS_Briggs_NGHP_17A in my data) are shown as 0. This is untrue according to QIIME2view outputs and core-metrics-results (generated following the moving pictures tutorial) and a barchart created using the phyloseq object. Additionally, some stats (alpha diversity) are calculated accurately yet other stats (core members of the microbiome) return errors citing the otu_table.
Ideal solution(s): 1) the otu_table from the phyloseq object is accurate, 2) I debug stats issues, and/or 3) all of the above.
Please see the following code and files for reproducibility below:
library("qiime2R")
library("phyloseq")

pyro_physeq<-qza_to_phyloseq(
  features="/path/pyro_sg.table.qza",
  tree="/path/pyro_sg.r.tree.qza",
  taxonomy = "/path/pyro_sg.taxonomy.qza",
  metadata = "/path/pyro_sg.met.short.txt")

#section out otu_table for ease
otu_table <- (otu_table(pyro_physeq))
otu_table #all SRRs are the study "AS_Briggs_NGHP_17A" except for the first two, SRR1539064 and SRR1795474

I am not expecting the majority of these values to be 0's. I can create barcharts with the same phyloseq object, which use the otu_table as the input.
plot_bar(pyro_physeq, fill="Class")

Alternatively, I've tried to section out the study "AS_Briggs_NGHP_17A" to visualize alone, and the barchart is how I'd expect. Similarly, I tried a quick shannon diversity analysis which was also calculated (and returned believable values).
phy_briggs <- subset_samples(pyro_physeq, name == "AS_Briggs_NGHP_17A")
plot_bar(phy_briggs, fill="Class")
div.sh <- data.frame(diversity(otu_ord, index = "shannon"))

HOWEVER, if I extract that otu_table it's all 0's. If I attempt to calculate the core microbiome using 'core_members' (R package microbiome), the 0's for all ASVs return errors.
pseq.core <- core_members(phy_briggs,prevalence = 50/100,detection = 1/100)
pseq.core
# returns 'character(0)', or to see the error:
pseq.core <- core(phy_briggs,prevalence = 50/100,detection = 1/100)
Error in validObject(.Object) : invalid class “otu_table” object:
 OTU abundance data must have non-zero dimensions.

This also occurs for the full phyloseq object (pyro_physeq). Thoughts?
Data
metadata: https://www.dropbox.com/scl/fi/wloiq8dze3znplmoqytzt/pyro_sg.met.short.txt?rlkey=bfsnrdqt1peh52xc2cqb3tgnd&dl=0
tree: https://www.dropbox.com/scl/fi/oliaj1qb9annxxptlmuix/pyro_sg.r.tree.qza?rlkey=ijmjeyt77tvk63cn1tovbcgfr&dl=0
feature table: https://www.dropbox.com/scl/fi/goav68y9osdhu6zoam2e0/pyro_sg.table.qza?rlkey=vzmbksfcpf1wvfyqy1frozxcs&dl=0
taxonomy: https://www.dropbox.com/scl/fi/3nmg3aobl8ck04uvm827u/pyro_sg.taxonomy.qza?rlkey=ujs4tjtzih5i6rhzj9xl4icri&dl=0

Thanks for the detail. Is this a 16S tree? Is the underlying data set 16S metagenomics or whole genome metagenomics?

Hello, it's all 16S rRNA sequencing using primer-tagged pyrosequencing.

I do understand microbial metagenomics in detail. I would just print out the tree and see what is there. ggtree should do that, i.e. for pyro_physeq. Maybe just post the treefile and see if it looks reasonable. I appreciate the detail of the question, but to fully resolve this appears to be a lot of work.


text:
Question about umap using different numbers of pca components as initialization
I am new to the scRNA-seq field and I have been doing some experiments of visualization of UMAP using different numbers of PCA components for initialization. The process involves projecting scRNA-seq data (count matrix) onto various numbers of PCA components, followed by non-linear dimension reduction using UMAP. It seems that the UMAP results are exactly the same across different initializations, regardless of the number of PCA components used. Are there any theoretical reasons for this? I am using Dynamo and their function dyn.tl.reduceDimension(adata, n_pca_components=30) for the experiments.
The results can be reproduced using the following code. Here I used the build-in zebrafish dataset with n_obs × n_vars = 4181 × 16940. I am just trying to see how the PCA initialization affects the Umap results.
import warnings
warnings.filterwarnings('ignore')

import dynamo as dyn

dyn.configuration.set_figure_params('dynamo', background='white')

adata = dyn.sample_data.zebrafish()

preprocessor = dyn.pp.Preprocessor(cell_cycle_score_enable=True)
preprocessor.preprocess_adata(adata, recipe='monocle')

dyn.pl.variance_explained(adata)

dyn.tl.reduceDimension(adata, n_pca_components=1, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

dyn.tl.reduceDimension(adata, n_pca_components=30, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

Here is the picture that shows the variance carried by PCA:

And the final umap results using n_pca_components= 1 and n_pca_components= 30:

Thank you for your reply @M__. I have tried actually n_pca_components= 1, 2, 5, 10, 30, and all umap results are the same. More information has been added.

Thanks for your edits and new information.

It seems you are performing PCA -> UMAP. The results you are reporting actually look cool. It's a good answer and there's nothing to worry about. It also shows good analytics, 'cause you're questioning the robustness of your data.
The limiting factor is the variance carried by PCA - which needs to be calculated (it may be low, e.g. <50%). Well I mean you've calculated it via the PCs - but its not declared, its better to know how much it is, i.e. how much variance each PC holds. The UMAP is doing nothing more than representing the established variance identified by PCA.
If you change the number of components which UMAP is assessing then this will result in changes in the clusters observed, usually.
However if you are changing the no. PCs -> UMAP and there's little effect then the PCs are not introducing new information to cut-up the data set. Thats a fairly robust answer: it's attained "stationarity". Very desirable output, just make sure there isn't a bug in your code, but it demonstrates there's no hidden or partial information that might give a different result. So the clusters you are seeing are a good representation of the data (thats good).
Its not a complete analysis - as hinted - but I'd be happy with that result.

Overall, there's very heavy dispersion of the variance between PCs. I've never had to teach this stuff, beyond a few line explanations. An educator experienced with this theory would give a more succinct explanation. However, the discussion below arrives at a common understanding between the OP and I. If there's demand I'll summarise the discussion as a singular paragraph: the final comments should be the introduction, i.e. what PCA is doing at basic level.
Comments (discussion)
@Zach even when using n_pca_components=1 and n_pca_components=30, the results remain unchanged. The first PCA component should only capture a part of the overall information. Is this correct?
Response Number 0 (on your graph) is the first principle component (16%). PC = 1 is the second PC (5%). Thus PCA = 1 means 21% of the variance. This far exceeds the PC 3 at 2%. Thats why it works, so all information is PC1 and PC2, nothing beyond that discriminates against the data.
@Zach So, n_pca_components = 3 could capture almost all the information and adding more components doesn't really make any difference to the results. But if we only consider n_pca_components = 1, why there are still no differences compared to n_pca_components = 30?
Response its easy .. because if you add all the PCs from PC 4 + ... + PC 30 it will be <5% of the total variance. PC1 + PC2 + PC3 = 24% of the variance, so there is nothing to see beyond PC3, when the maximum is PC 30.
@Zach Why there are no differences between n_pca_components = 1 (using only PC1) and n_pca_components = 30 (using PC1+PC2+...+PC30)?
Response
PC1 + PC2 = 24% and PC1 ... + PC30 <30%. There is no information in the variance for between group clusters for PC3 to PC30, there is likely to be within group differences but you are not measuring those.
@Zach there is little difference between PC1 + PC2 = 21% and PC1+...+PC30. I mean if we only consider PC1 = 16% and PC1+...+PC30 > 25%, do they have any differences between PC1 and PC1+...+PC30?
Response PCA works on individual elements it, it doesn't work on groups. The analysis of groups is post-hoc (usually via colour coding). PC=1 in your UMAP calculation is PC1 + PC2 (PC1 is automatic, its position 0 in your graph). . You are only looking at between group variance and in your data that is resolved essentially for PC1 + PC2 = 21%. The within group positioning will change so, yes there will be difference against PC1+PC2 vs PC1 + .... + PC 30 but we can't see them. Again PCA is individual elements, it's not using group information. So if you did sub-group analysis you might see changes PC1+PC2 vs PC1 + ... PC 30. So to put it simply your group delineation is too big to spot in any differences (but that only matters if it is biologically relevant).
@Zach Actually what I mean is PC0 = 16% and PC0+...+PC30 > 25%. But your explanation makes sense to me and I appreciate your response.

Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on [meta], or in [chat]. Comments continuing discussion may be removed.

It is a bit odd that different numbers of principal components are giving identical results, especially given that even the same initialisation settings with different UMAP epoch numbers will give different results.
I'm not familiar with the use of python for single cell data analysis, but it looks like the adata information might be being modified when running reduceDimension. If this is the case, then it's possible that the first PCA reduction is reducing the dataset as far as it will go, and the second PCA reduction doesn't do anything. This could be checked by changing the order of operations (from a fresh run) and seeing if the output is changed:
dyn.tl.reduceDimension(adata, n_pca_components=30, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

dyn.tl.reduceDimension(adata, n_pca_components=1, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

It's also possible that the umap command is doing its own reduction and ignoring the result of dyn.tl.reduceDimension. Maybe the UMAP plotting function doesn't actually do a UMAP reduction as well, and the UMAP result (or initial condition) is calculated earlier (e.g. in preprocessor.preprocess_adata).
It's also possible, as @M__ suggests, that you've reached a stationary result after however many UMAP iterations are set as default. But... I find that unlikely given how variable UMAP has been for me. I'd want to do a lot more exploration before settling on that as an explanation.


text:
What is happening in the Zuker-Algorithm
I am trying to understand the Zuker Algorithm.
What is happening during initialisation? I know in Nussinov the diagonal is filled with 0s. Are the diagonals filled with infinite energy during Zuker´s initialisation?
Also, what is happening during recursion?
I know unlike Nussinov, Zuker pays attention to stacking energies and loop structures, not individual basepairs. But how? How does the algorithm calculate the energies for stacking and loop structures at position 1 if there is no info about its neighbors yet?
I understand there are 2 matrices - V and W. Is it correct that V is part of W? And is it correct that V keeps information about stack interactions and how much energy a stack of base pairs contributes - like keeping score on a multiplier?
Is there a code or example (I would prefer R), which demonstrates how the algorithm goes over a sequence? Like a demonstration. I thought of one myself, but I am not sure if its correct:
sequence <- c(1, 2, 3, 4, 5,6,7,8,9,10,11,12,13,14,15)

n <- length(sequence)

for (l in 2:n) {
   for (i in 1:(n - l + 1)) {
      j <- i + l - 1
      cat("i =", i, "j =", j, "\n")
   }
}

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 am using RNAFold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). I am suppose to write a shor text after my work and I am having diffuculties understanding certain espects of Zuker`s Algorithm. Nussinov is easy to understand and visualize on paper. And i was wondering how would implement Zuker in R or Python for example or how one would fill a dotplot using Zuker.

May i ask.. Whats the graph theory code in Python. I found some Python code on Github which is suppose to simulate Zuker, but if its not the correct Zuker Algorithm, then I would like to know.

Graph theory code for RNA structure is a separate question. This question is about Zuker's code - I don't know it - I thought Zuker used graph theory (which I know) but used matrix algebra instead. It's easy to get hold of the code just download the src for MFold and especially RNAfold. RNAfold is in C.


text:
How to access yeast-cyc in a programmatic way to extract the GO terms?
I am a newbie in Bioinformatics, I want to extract the GO terms from yeast-cyc,
which is items in the Go terms tab of this page https://yeast.biocyc.org/gene?orgid=YEAST&id=G3O-29622#tab=GO , not one page but one database like Saccharomyces cerevisiae S288c.
Is there an API or python lib can help us finish this task? I had found biocyc in pypi, but I did not found Go terms was involved in this API.
If the API is not convinient, can I download some dataset in which my required items is include?

Not sure if I've misunderstood, but if you want the GO terms for yeast you can get them from https://www.uniprot.org/uniprotkb?query=(taxonomy_id:559292). Using the UniProt system you can then generate an API endpoint URL for your query, e.g. here are the GO terms for all yeast genes: https://rest.uniprot.org/uniprotkb/stream?fields=accession%2Cid%2Cprotein_name%2Cgene_names%2Cgo&format=tsv&query=%28%28taxonomy_id%3A559292%29%29 (use wget or cURL to pull the data); does that help?

Thanks, it is the data I need.

Great! I'll turn my comment into an answer :)

You can get the GO terms for yeast from UniProt, e.g. https://uniprot.org/uniprotkb?query=(taxonomy_id:559292).
Using the UniProt website you can generate an API endpoint URL for your query, e.g. here are the GO terms for all yeast genes: https://rest.uniprot.org/uniprotkb/stream?fields=accession%2Cid%2Cprotein_name%2Cgene_names%2Cgo&format=tsv&query=%28%28taxonomy_id%3A559292%29%29
You can download this data using wget or cURL, e.g.
wget -O yeast_go_terms.tsv "https://rest.uniprot.org/uniprotkb/stream?fields=accession%2Cid%2Cprotein_name%2Cgene_names%2Cgo&format=tsv&query=%28%28taxonomy_id%3A559292%29%29"


text:
Finding mutations in glycosylation sites
We are going to look at the likely N- and O- glycosylation sites within MUC16 (Q8WXI7 · MUC16_HUMAN)§ by in house long-read DNA sequencing data (PacBio).
Counting the number of tandem repeats is a task which could be done with bioinformatic tools like ad-codeVNTRas shown in the below figure. Another task is identifying mutations within these repeats. My question is: after predicting glycosylation sites using tools like PA2DE within these tandem repeats in MUC16 gene, how could I identify any point mutations in these glycosylation sites?
I googled a lot to get some idea; Shall we do germline analysis in MUC16 glycosylation predicted sites? Is this what we should use to look at glycosylation sites that are ‘rejected’ to understand the sites that can’t be glycosylated (mutated)?
Does anyone have such an experience how to identify mutation in glycosylation sites in MUC16 which known to have variable number tandem repeat (VNTR)?

If you write the exact rules, it's not difficult to write a regex function to find these. How many repeats? What does tandem mean? Is it without spaces, VNTRVNTRVNTR?

Thank you, counting the number of repeat is another task which could be done with bioinformatic tools like ad-codeVNTR. Another task is identifying mutations within these repeats. My question is after predicting glycosylation sites, how could I identify any point mutations in these glycosylation sites.

Please [edit] your question and clarify what you need and what you are working with:  What are your input data? Protein? DNA? Have you performed your own sequencing? Are you looking through published sequences? How do you define mutation, what is your reference? Doesn't code-adVNTR do exactly what you want? The image you show suggests it does, so isn't it already what you need?

Thank you I tried to edit my question. That would DNA long read sequencing done on human cancer cell lines


text:
What options are there for visualising DNA sequences?
I typically see several ways of visualising DNA sequences and I would like to know when (best) to use them?
For example, I sometimes see the letter codes GATC, sometimes I see colours, sometimes I see a representation of the double-helix. I think I saw a circular representation.
Additionally, what other options exist for visualising DNA and when should they be used? E.g. for papers, analysis in software, etc.
What do people tend to use when presenting DNA sequences?

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

Hi @ScottMcKenzie. Welcome! Hope you stick around. There are loads of ways to visual sequences and loads of different sequence formats, e.g. alignments, assemblies (NGS stuff) with annotation, nucleotide/allele frequency plots (big letter stuff) ... there are genome and vector (non-mathematical) representations. Could you kindly narrow the question to your specific aim/area please? Maybe tell us your project, type of DNA you want to represent? Otherwise this question will become a thesis of answers.

Yes, the best way will depend entirely on what you are trying to present. Why are these sequences relevant? What are you trying to say about them?

DNA sequences have a massive range in length, so the approach of visualisation is dependent on the purpose of that visualisation.
I have been experimenting with visualising the repetitive structure of DNA in a way that represents repeats at both short range (tens of bases) and long range (millions of bases, or longer). You can see my attempts at that here. An associated presentation about DNA visualisation (with emphasis on repetitive sequences) can be found here.

This is brilliant, thanks. PS I'm in Auckland! Small world... :)

If one of the answers here solved your issue, @ScottMcKenzie, please take a moment and accept it by clicking on the checkmark on the left. That is the best way to express your thanks on the Stack Exchange sites.


text:
How to analyse qualitatively the penetration ability of particles in spheroids using fluorescent z-stacks?
0
I have to establish the penetration profile of particles in tumour spheroids.
For this I have spheroids composed of cells which were exposed to particles, both of which are fluorescently labeled. The spheroids were then imaged using a confocal microscope by taking z-stacks of all spheroids. For analysis, I determined the intensity profile of each focal plane using an ImageJ plugin, namely "distribution_shapes".
What should be done next?
I currently did not find any step-by-step guide for this analysis for manual assessment and any type of help would be greatly appreciated.

This looks like an interesting image analysis question, which is great because we don't get as many of those here. Could you be more specific about what you mean by "penetration profile"? Can you point to any publications that use a similar analysis? I find one possibly related paper here but they don't go into any detail unfortunately.

The next step depends heavily on what your question is. Some general ideas, not specific to tumor particle penetration.

Your current visualization seems OK to me, if you want to make some features more obvious you can try density estimation.

If you have a generative model (for example you believe this is a mixture of a "shallow" and "deep" Gaussian distribution, or you have reason to believe the particles follow an exponential distribution), you can fit your data to this distribution and estimate its parameters.

If you want to compare conditions without too strong assumptions on the distributions, you could compute the cumulative distribution of the particles. There are statistical approaches to compare such curves.

To compare conditions, you might even be able to come up with a simpler metric, for example "total fluorescence in the core of the spheroid" normalized by "total fluorescence in its periphery".

Ultimately it all depends on your research question (are there treatments to compare to each others, biological replicates whose consistency you want to check, known structure within the tumor spheroid, a time course, ...). My suggestions above would be easier to implement in R/Matlab/Python than pure ImageJ.

I want to compute the cumulative distribution of the particles. But I don't know how to do this with a z-stack as the distribution across z-direction will not be the same as in x and y-directions. I also don't have any knowledge of programming, that's why I was wondering if anyone knows a guide or course that can help in this regard. I did find a matlab script, but I get the same errors and I asked several people to help me with that and we didn't manage to solve it.

I am not sure if this is exactly what you're looking for, but I was able to find this ImageJ script that estimates a penetration profile in tumor spheroids.
I am not an ImageJ person but I suspect that with some minimal changes this could be used to answer your question.

This is actually what I was using to determine the intensity profile of my particles. However, I am not sure what to do next, because I have the intensity profile of a spheroid over z-axis and don't know how to put the information together or to visually represent the quantification of the penetration profile.

@Timi thanks for clarification. Can you post what your data look like (just a snippet is probably enough)? That might be something we can help with.

I added a screenshot of how to profile intensities look across a z-stack (1st image) and how the average profile intensity looks like in the last 5 focal planes in the z-stack. This is what we are aiming to get.
The papers that I found that do something similar are: https://www.sciencedirect.com/science/article/pii/S0168365921006672
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7824314/
, but they also don't go into details, and the script in Matlab is not working for me.

@Timi Can you show the original raw data (text, not plots) output of the imagej protocol, rather than plotted summaries? As in, what is the data that you used to make those plots. It seems that you want to make better or more informative plots similar to some that you are seeing in those papers (which figures? maybe include them in the question). Just a subset of the data is fine, to show what the format is, or even faked data that is in the same format if you don't feel comfortable posting raw data. I don't see the matlab script you are talking about, a direct link might help.

I don't know how to attach a table to this post. Would it be possible to send you an Excel sheet in an email or private message?

@Timi What I'd suggest doing is using table markdown or similar to embed a table in your question. There are web converters for converting spreadsheet to markdown. Alternately you could use a google drive link, pastebin, dropbox, box, etc. In each case I'd suggest embedding the link directly in your question so that anyone who needs to can look at it.


text:
How should I use dorado basecaller to calculate translocation time? Can I change models' config files?
I am trying to use dorado basecaller to process ONT data. I need to make the move table and I am doing it with --emit-move.
dorado basecaller dna_r9.4.1_e8_sup@v3.3/ /pod5_files/ --recursive --emit-moves --modified-bases 5mCG_5hmCG -v > /bams/mod_stride1.bam

The main purpose of the move table for me is to calculate the number of ONT device measurements before basecaller emits one base (calculate translocation time per base). But the stride value in config.toml file is fixed to 5 (stride=5). I need to calculate the instantaneous translocation time per base (in seconds), because I am looking to find a measure similar to IPD values in PacBio, when we can get the time for sequencing each base.
I was wondering to change the stride value in config.toml file in the model folder, from stride=5 to stride=1, to directly infer the number of measurements without using Tombo for correcting and re-squiggling.
Unfortunately, even though this command works perfectly fine when the stride=5, when I set stride=1 it always gives me the following error: "CUDA out of memory", even when I triple the size of RAM for the job.
My questions are:

Does setting stride=1 make sense for this purpose? Or do I have to do something else? (I know we have some K-mer models too)
Is this even possible to change the stride value on config.toml file?
How should I figure out how much RAM I need for the job?

In pod5 files the number of samples prior to the adapter/sequence should be encoded in the metadata, and can be looked at using ONT's pod5 inspect tool:
$ pod5 inspect reads pod5_file.pod5

read_id,channel,well,pore_type,read_number,start_sample,end_reason,median_before,calibration_offset,calibration_scale,sample_count,byte_count,signal_compression_ratio
00445e58-3c58-4050-bacf-3411bb716cc3,908,1,not_set,100776,374223800,signal_positive,205.3,-240.0,0.1,65582,58623,0.447
00520473-4d3d-486b-86b5-f031c59f6591,220,1,not_set,7936,16135986,signal_positive,192.0,-233.0,0.1,167769,146495,0.437
...

The average translocation speed can be worked out by comparing the sample_count to the number of called bases.
More details here.
Unfortunately, due to the way the basecaller models work, I don't think there's any easy way to create a 1-sample lookup table for precise instantaneous translocation calculations.
If it is just an out of memory error, then it should work (but very slowly) if you set the basecaller batch size to 1 (i.e. -b 1). If even that doesn't work, something else is going on and it's unlikely to be a memory problem.


text:
Ordering average of gene expression in heat map
I want to plot a heat map in header I have the average of gene expression validated in lab for each sample
dput(heat[1:5,1:5])

structure(c(3.05574677859363, -0.451189292623471, -0.398066832453666,
-0.456927931251448, -0.393418432279742, -0.270234667496526, -0.451189292623471,
-0.398066832453666, -0.456927931251448, -0.393418432279742, -0.270234667496526,
-0.451189292623471, -0.398066832453666, -0.456927931251448, -0.393418432279742,
-0.270234667496526, -0.451189292623471, -0.398066832453666, 2.0338091136055,
-0.393418432279742, -0.270234667496526, 1.76078923919148, 1.96896026174603,
1.97960246116808, -0.393418432279742), .Dim = c(5L, 5L), .Dimnames = list(
    c("hsa-miR-548b-5p", "hsa-miR-181c-5p", "hsa-miR-181d-5p",
    "hsa-miR-125a-3p", "hsa-miR-4286"), c("30566-001", "30566-002",
    "30566-003", "30566-004", "30566-005")))

dput(m[1:5,c(2,7)])

structure(list(HealthStatus = c("Healthy", "Healthy", "Healthy",
"Healthy", "Healthy"), Average = c(108.915677, 103.148873, 40.3359607,
88.4329528, 87.7756577)), row.names = c("30566-001", "30566-002",
"30566-003", "30566-004", "30566-005"), class = "data.frame")

pheatmap(heat, cluster_rows=T, show_rownames=T,
         cluster_cols=FALSE, annotation_col=m[,c(2,7)],
         annotation_colors=annot_colors)

How I can order average in a descending way?

If you order your matrix/data.frame you can sort your columns by decreasing "Average", e.g.
library(pheatmap)

heat <- structure(c(3.05574677859363, -0.451189292623471, -0.398066832453666,
                    -0.456927931251448, -0.393418432279742, -0.270234667496526, -0.451189292623471,
                    -0.398066832453666, -0.456927931251448, -0.393418432279742, -0.270234667496526,
                    -0.451189292623471, -0.398066832453666, -0.456927931251448, -0.393418432279742,
                    -0.270234667496526, -0.451189292623471, -0.398066832453666, 2.0338091136055,
                    -0.393418432279742, -0.270234667496526, 1.76078923919148, 1.96896026174603,
                    1.97960246116808, -0.393418432279742), .Dim = c(5L, 5L), .Dimnames = list(
                      c("hsa-miR-548b-5p", "hsa-miR-181c-5p", "hsa-miR-181d-5p",
                        "hsa-miR-125a-3p", "hsa-miR-4286"), c("30566-001", "30566-002",
                                                              "30566-003", "30566-004", "30566-005")))

m <- structure(list(HealthStatus = c("Healthy", "Healthy", "Healthy",
                                     "Healthy", "Healthy"), Average = c(108.915677, 103.148873, 40.3359607,
                                                                        88.4329528, 87.7756577)), row.names = c("30566-001", "30566-002",
                                                                                "30566-003", "30566-004", "30566-005"), class = "data.frame")

pheatmap(heat[,order(m[,"Average"], decreasing = TRUE)],
         cluster_rows=TRUE, show_rownames=T,
         cluster_cols=FALSE, annotation_col=m)

Created on 2023-09-06 with reprex v2.0.2


text:
How to analyze GC content of fasta file?
I have the assembly file of a certain lizard species, I wanted to recreate Fig. 2 of the following paper, which compares the GC content distribution over 5kb windows, how would I go about doing this?
https://bmcecolevol.biomedcentral.com/articles/10.1186/s12862-022-02086-7#Sec10

I would suggest using bedtools makewindows and bedtools nuc to compute the GC content for genomic windows as shown in Fig 2 of the linked paper.
See related question here.
This will slice the genome into 5Kbp (or whatever size specified to makewindows) pieces and compute the GC content of each piece, which can be plotted as a distribution as in the figure.
Something like this would work for the 5000 bp example:
# assemble a bedtools compatible genome file:
# see https://www.biostars.org/p/70795/ for hints how to do this.
bedtools makewindows -w 5000 -g my_genome.genome > windows.bed # you just made this .genome file manually
bedtools nuc -bed windows.bed -fi my_genome.fasta > windows_gc.bed

The bed output file will contain the GC content of each window in the 7_pct_gc column (according to my understanding of the output file).

Please note that you do not have to make ".genome" file manually if you have samtools. You could simply generate an index file by running samtools faidx genome.fasta, and use the output file with extension .fasta.fai as input in your makewindwos command.

tail -n +2 windows_gc.bed | cut -f 5 | awk '{printf("%d\n",$1*100)}' | sort | uniq -c | sed 's/^ *//'| sed 's/ /,/'>output.csv. This should print a CSV file which you can easily plot.


text:
Why does the subset of rows that match a condition added by the rows that don't match a condition not add up to the total number of rows?
I need to match the external gene names of a table I have called (which has 1,000,000 entries) with their respective ensemble gene IDs.
Here's a small subset of df:

V1
external_gene_name
V3

M0085
A1BG
1

M0085
SNX8
1

and here's a small subset of ensembl_ids:

ensembl_gene_id
external_gene_name

ENSG00000288263
GRK1

ENSG00000291428
BRSK2

I also included a line to see the number of entries in external_gene_name which do not have corresponding ensemble names. Here is my code:
merged_df <- merge(df, ensembl_ids, by = "external_gene_name")
elements_not_in <- df$external_gene_name[!(df$external_gene_name %in% ensembl_ids$external_gene_name)]

But what's confusing me is that when I add the number of rows in merged_df with that of elements_not_in, the total exceeds 1,000,000. Shouldn't the total be 1,000,000? I'm wondering if there's a flaw in my code.

From a brief glance, the gap is not in your code but in your expectations. You're operating based on an assumption of 1:1 mapping between ENSG IDs and HGNC symbols (your external_gene_name attribute). In reality, multiple ENSG IDs will map to the same HGNC symbol so the merge will produce multiple rows per external_gene_name in the merge product. This is why you're seeing a total greater than the sum of number of unique entries when you compare the merge product + the anti-join result (i.e., the elements_not_in dataset) to the initial individual datasets.


text:
How to perform PCA on proteomic data set
I'm trying to perform Principal Component Analysis using R on a proteomics dataset.
As the dataset contains a lot missing values I tried different approaches.

I ran PCA using prcomp() and autoplot() after deleting all samples with missing values. This worked out nicely, but the plots are not representative as I had to delete more than 70% of the samples.

I imputed the missing values and ran a PCA (again prcomp() and autoplot()) which worked out, but provided me with a PCA plot showing a strong horseshoe effect (which seems off as the PCA I performed on the complete samples of my data doesn't show such a distribution).

I then came across this article (https://www.sciencedirect.com/science/article/pii/S1574954121000261) which introduces a modified PCA, InDaPCA (PCA of Incomplete Data). I transformed the code provided by the authors according to my data but again got errors/warnings because of missing values.

This is the structure of my dataframe, the code I used to perform InDaPCA as well as the resulting error:
num [1:261, 1:1132] 13.18 13.18 14.14 8.32 9.11 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:1132] "protein.A0A075B6H9" "protein.A0A075B6I0" "protein.A0A075B6I1" "protein.A0A075B6I9" ...

InDaPCA <- function(protein_data_filtered) {

  #scaling
  X <- scale(protein_data, center = TRUE, scale = TRUE)

  #correlation
  C <- cor(X, use = "pairwise.complete.obs")

  #Eigenvalue
  Eigenvalues <- eigen(C)$values
  Eigenvalues.pos <- Eigenvalues[Eigenvalues > 0]
  Eigenvalues.pos.as.percent <- 100 * Eigenvalues.pos / sum(Eigenvalues.pos)

  #Eigenvectors
  V <- eigen(C)$vectors

  #Principal components
  X2 <- X
  X2[is.na(X2)] <- 0
  PC <- as.matrix(X2) %*% V

  #object.standardized
  PCstand1 <- PC[, Eigenvalues > 0] / sqrt(Eigenvalues.pos)[col(PC[, Eigenvalues > 0])]
  PCstand2 <- PCstand1 / sqrt(nrow(PC) - 1)

  #loadings
  loadings <- cor(X, PC, use = "pairwise.complete.obs")

  #arrows for biplot
  arrows <- cor(X, PC, use = "pairwise.complete.obs") * sqrt(nrow(X) - 1)

  #output
  PCA <- list()
  PCA$Correlation.matrix <- C
  PCA$Eigenvalues <- Eigenvalues
  PCA$Positive.Eigenvalues <- Eigenvalues.pos
  PCA$Positive.Eigenvalues.as.percent <- Eigenvalues.pos.as.percent
  PCA$Square.root.of.eigenvalues <- sqrt(Eigenvalues.pos)
  PCA$Eigenvectors <- V
  PCA$Component.scores <- PC
  PCA$Variable.scores <- loadings
  PCA$Biplot.objects <- PCstand2
  PCA$Biplot.variables <- arrows

  return(PCA)
}

pca_result <- InDaPCA(protein_data_filtered)

Error in eigen(C) : unendliche oder fehlende Werte in 'x'
Called from: eigen(C)

I would highly appreciate any suggestions on how to perform PCA of datasets with many NA values and/or how to handle the horseshoe effect.
Thanks a lot in advance!

Not a direct answer to your problem but how about using a package tailored for proteomics datasets? For example msqrob2 workflow uses limma's plotMDS() after normalizing the data.


text:
Error output in nextflow pipeline using fasterq-dump
I have one file with multiple SRA accession called Dataframe_with_accession.txt, for this example I just put one SRA in that file called :SRR8933535
And the idea is to create a nextflow pipeline to download the sra files in the Dataframe_with_accession.txt, zip the files  using pigz and remove the previous fastq file
To do that I used this code :
params.inputFile = "/crex/proj/Dataframe_with_accession.txt"
params.outputDir = "/crex/proj/Output"

process DOWNLOAD_FASTQ {

    publishDir params.outputDir, mode: 'symlink'

    input:
    each (accession)

    output:
    tuple path("params.outputDir/${accession}_1.fastq.gz")
    path("params.outputDir/${accession}_2.fastq.gz")

    script:
    """
    fasterq-dump --threads 12 --outdir ${params.outputDir}  ${accession}
    pigz -p12 ${params.outputDir}/${accession}_1.fastq
    pigz -p12 ${params.outputDir}/${accession}_2.fastq
    """
}

workflow {
    // Read accessions from the input file
    accList = file(params.inputFile).readLines()

    DOWNLOAD_FASTQ(accList)
}

Then I run this nextflow file.
But I got the following error message :
[01/298f4b] process > DOWNLOAD_FASTQ (1) [100%] 1 of 1, failed: 1 ✘
Error executing process > 'DOWNLOAD_FASTQ (1)'

Caused by:
  Missing output file(s) `SRR8933535_1.fastq` expected by process `DOWNLOAD_FASTQ (1)`

Command executed:

  fasterq-dump --threads 12 --outdir /crex/proj/Output  SRR8933535

Command exit status:
  0

Command output:
  (empty)

Command error:
  spots read      : 8,389,148
  reads read      : 16,778,296
  reads written   : 16,778,296

Work dir:
  /crex/proj/Output/work/01/298f4be6dde2ee1856df31db169bec

Tip: you can try to figure out what's wrong by changing to the process work dir and showing the script file named `.command.sh`

And I do not undestrand what is going wrong here ? It says
Command output:
  (empty)

But the file is created

Are you aware of the fromSRA Nextflow channel factory? https://www.nextflow.io/docs/edge/channel.html#fromsra

Each nextflow process looks in its working directory for the files listed in output: . In your script: you already direct the output to its eventual destination. After running the script nextflow tries to find the output files and stops, because it can't find them. Solution: keep the output files in your workdir, nextflow will detect them and use the publishDir directive to copy  them to their eventual destination.

ok perfect thanks


text:
skip splits in character string
I need to generate a list of sequences from a database, in which a split pattern is skipped up to n-times. As I only know R, I would appreciate an R solution.
Example:
character string:
"GRGKA"
I could use a simple
strsplit()
function but this would result in a "perfect" split pattern.
What I am looking for is all combinations that are generated when n number of splits are skipped.
I specifically want to split the string after K or R (so that K and R are still included in the reported split sequence)
so a normal split would yield:
"GR" "GK" "A"
however with one "skipped" split it would look like this:
"GRGK" "GKA"
and with two skipped splits:
"GRGKA"
I was not able to find something online. Sorry if I missed something or did not provide enough info.
Here is the code that I started out with and which provides a normal split without skipping K or R
`"GRGKA" %>% strsplit(., "(?<=K|R)", perl = TRUE)`

[[1]]
[1] "GR" "GK" "A"