text:
Good / recommended way to archive fastq and bam files?
We have a lot of Illumina sequenced exome data. Currently we are using spring for its great lossless compression, but we are looking if there is anything better (and most preferably opensource) which can let us compress our fastq files.
We also want to compress or reduce file size of BAM, but using CRAM in lossless mode didn't yield good results. The space savings were decent, roughly ~10% (IMO was expecting more), but the biggest hit came when using it with IGV. We have set up a http-server and stream bams to IGV (Helped us immensely as people wouldn't have to download and visualize bam files). With Cram, the loading times were significantly high, requiring users to commit more memory for IGV. Sometimes we would have to wait around a minute for reads to show up. With Bam it is near instantaneous.
Can't use lossy compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).
What are some good ways to compress fastq and bam files so we can maximize our storage?

EDIT: I am rewriting the answer in response to updates to the original question.
TL;DR: use CRAM
Background 1: quality binning and FASTQ compression
In the old days, base callers outputted base quality at full resolution – you could see quality from Q2 to Q40 in full range. As a result, quality strings were like semi-random strings and very difficult to compress. Later people gradually realized that keeping base quality in low resolution wouldn't affect downstream analysis. The Illumina basecaller started to output quality in 8 distinct values and later changed that to 4 bins. This change greatly simplified quality string and made them compressed better. For example, in old days, a 30X BAM would take ~100 GB. With quality binning, it would only take ~60 GB.
Background 2: GATK Base quality recalibration
In early 2010s, Illumina base quality was not calibrated well. GATK people introduced BQSR to correct that and observed noticeable improvement in SNP accuracy. Nonetheless, with improved Illumina base caller, their base quality became more accurate. Meanwhile, the world moved to 30X deep sequencing. The depth overwhelms slight inaccuracy in quality. I would say around 2015, BQSR was already unnecessary for data produced at the time.
Does it hurt to apply BQSR? Yes. First, BQSR introduces subtle biases towards the reference and towards known SNPs. Second, BQSR distorts the data. At least for some datasets, I observed that SNP accuracy dropped with variant quality after BQSR; I didn't observe this with raw quality. Third, BQSR is slow. Fourth, for new sequencers producing data at higher quality, BQSR is likely to decrease data quality. Last, related to the question, BQSR added another semi-random quality string and made compression even harder. Nowadays, running BQSR is a waste of resource for worse results. The official GATK best practice no longer uses BQSR according to their WDL file.
CRAM
These days a 30X human CRAM only takes ~15 GB (see this file). This is a huge contrast to ~100 GB BAM in early 2010s. OP only saw ~10% saving probably due to a) BQSR and/or b) old data with full quality resolution.
On encoding/decoding speed, CRAM was much slower than BAM. Not any more. The latest htslib implementation of CRAM is faster than BAM on encoding and only slightly slower on decoding. The poor performance of IGV on CRAM could be that the Java CRAM decoder is not as optimized.
It is true that CRAM is not as widely supported as BAM. However, all the other alternatives are much worse. Petagene said they had IGV-PG (PDF) for their format. That is not official IGV and I couldn't find more recent update beyond the 2019 press release. I don't see other viable options.
Note that the common practice is to keep all raw reads, mapped or not, in BAM or CRAM such that you can get raw reads back later. BAM/CRAM additionally keeps metadata like read group, sample name, run information etc and is actually more popular than FASTQ in large sequencing centers. Also note that you don't need to sort CRAM by coordinate. Unsorted CRAM is only a little larger than sorted CRAM.
CRAM and its competitors
The core CRAM developer, James Bonfield, is one of the most knowledgeable researchers on compression (and one of the best C programmers) in this field. He has done a lot of compression evaluation over the years. The conclusion is that on a fair benchmark, CRAM is comparable to the best tools so far in terms of compression ratio.
Petagene could compress better in the plot @terdon showed mostly because it has a special treatment of the OQ tag generated by GATK BQSR. It is a typical trick marketing people use to make their methods look better. With BQSR phased out, this plot is no longer relevant.
On commercial software
In general, I welcome commercial tools and think they are invaluable to users. I also have huge respect to DRAGEN developers. However, on FASTQ storage, I would strongly recommend against closed-source compressors. If those tools go under, you may lose your data. Not worth it.

Thanks for your input, but mainly for the link. WDL seemed too complicated for me. We also noticed that BQSR was not recommended by nf-core people as well, but none provided a concrete reason beyond "its not useful". As I've also commented on @gringer post, CRAM didn't work well when using it with IGV, as well as some downstream tools don't work with CRAM.

@KarthikNair It is interesting that you are criticizing CRAM for slowness in IGV, while you are compressing fastq with spring and you can't even load other alignment formats into IGV. What are you expecting? Also, you don't need to understand WDL. Just search/grep "BaseRecalibrator" or "BQSR" to see if it gets called.

@KarthikNair sadly, the GATK best practices, which were updated very recently, still recommend BQSR so it is really surprising that their recommended pipeline doesn't use it. That said, in my place of work, we analyze around 2k human samples a month, and we have found that BQSR often causes us to miss variants. I don't have statistics for this yet so it's anecdotal, but there is strong evidence that BQSR not only isn't particularly helpful today, but can be actively harmful. And not only with small panels, this also happens with whole exome sequencing samples.

@terdon That page is referencing their older pipelines, which do use BQSR. The WDL file in my link is two years old. It is the latest I can find. IMO, BQSR has been hurting for 10 years. Well, better late than never.

you focus on BAM but what would you suggest for FASTQ if not DRAGEN ORA reference based approach? Arguably long term of storage of FASTQ is more important since BAMs can be regenerated

@Chris_Rands See new edits. I see the value in reference-free algorithms such as spring and dsrc. But if you already require a reference, it is better to keep CRAM than ORA.

For what it's worth, petagene's compatibility layer lets both IGV and JBrowse work out of the box with peta-compressed bam files. You don't even need to use their IGV implementation, vanilla IGV works fine as long as you are using petagene's trick with the symlinks I mentioned in my answer.

@Chris_Rands You can align raw reads with STAR or any local aligner and then create a CRAM. I guess that will give you a smaller file than BAM.

Hay all, thanks for all your comments! I will echo @terdon and encourage people to share their findings w.r.t. to BQSR not being very helpful. We have been using MarkDuplicatesSpark and BQSR pipeline. I'm planning on dropping the BQSR steps and using samtools markdup -s instead of MarkDuplicatesSpark (Mainly for optimization). And yes, perhaps the lackluster size reduction could have been due to BQSR in data. We also have old NGS data sequenced from 2016 onwards, so some samples might use unbinned quality scores.

I'd recommend reconsidering CRAM. It's an open standard that is most likely to still be accessible many years in the future. Regarding your point about lossy compression...

Can't use lossy [CRAM] compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).

By default, CRAM compression is not lossy.
CRAM does, however, have an optional archive settings mode (samtools view --output-fmt cram --output-fmt-option archive ...) which is a lossy compression, doing things like removing read names, removing additional accessory fields, and additional compression of quality scores. In all cases, the base sequence of the reads is preserved.
https://www.htslib.org/doc/samtools.html

You might want to consider PetaGene. This is a paid tool, and not open source, but it is something we use at work and it is really good! It can do lossless compression of gzip-compressed fastq files, bam files, and even cram files, and achieves really good size reduction (image from the PetaGene webpage):

In my own testing, using files of different sizes, I get an average size reduction of 62.1% (2.6 x smaller) for fastq.gz files (compared to the gzipped fastq file, not compared to the original, uncompressed fastq file) and an average reduction of 82.1% (5.6 x smaller) for bam files. I didn't test cram files.
If you have the budget for a paid service, these guys do good work! The greatest advantage, apart from the very significant size reduction, is that they also provide a system that allows you to use the compressed files natively. So you can feed the petagene-compressed bam file directly to bwa for example. They create special symlinks (foo.fastq.gz -> foo.fastq.petacompressed) and you can then use these symlinks and petagene handles the decompression on the fly.

The Wikipedia page on CRAM has a nice table with a few alternatives to CRAM which I reproduce below. I haven't tried any except PetaGene though so cannot comment on their performance:

Name
Key advantage over CRAM
Repo

Deez
None. An early approach to SAM compression, introduced 2014.
https://github.com/sfu-compbio/deez

Genozip
Significantly better compression, but slower random access
https://github.com/divonlan/genozip

GenomSys
Utilizes ISO-standard MPEG-G, but compression inferior to CRAM 3.1
Source code not publicly available

PetaSuite
Virtual filesystem interface
Source code not publicly available

We considered this software, but sadly we don't have the budget for this. And yes, on-the-fly decompression is a great feature to have, makes storing all data in this format extremely viable. Right now, we get decent compression with spring, but it needs lots of de/compression time, not to mention the CPU resources it guzzles during the process lol

@terdon genozip is quality work and is very impressive. Its full source code is also on github. If genozip goes under some day (I hope not!), you may create a fork and continue to use their format.

@user172818 ~terdon - can you share any further info on how PetaGene achieves compression and the critique mentioned of "special treatment of the OQ tag generated by GATK BQSR"? I know it's closed source but trying to understand conceptually how their compression works

Not I, @Chris_Rands, sorry. I'm sure user172818 will know far more, but I have no idea and wasn't even aware there was any special treatment before reading their answer. I know at least one person who used to be an engineer working for them is active on the site. If there's something that can be shared, maybe they will answer if you post a question.

An alternative to the good option of PetaGene that @terdon mentioned and also not open source is the DRAGEN ORA compression. This is loseless reference based compression and can be ~5 x smaller over fastq.gz, but it's only for fastq files and only those generated by Illumina machines (and will be easiest within the DRAGEN ecosystem).
Another point to consider is where you store the data. Even with 'standard' compression you might get more benefits moving archival data to a system like AWS S3 Glacier, depending on your current infastructure setup of course

This is exactly what I do. I archive FASTQ.GZ (with added content-md5 checksum) to S3. I also archive processed outputs to a different location that is less fault tolerant than S3 as those files can be recreated.

We did check this out, but we do not sequence the NGS data ourselves. We ask external labs to perform the sequencing and they share the fastq.gz files with us. What is the approximate pricing for this? I can't clearly identify where those details are!

The only free and open source tool I know that can help is zstd. Their github repository's README describes it as:

Zstandard, or zstd as short version, is a fast lossless compression algorithm, targeting real-time compression scenarios at zlib-level and better compression ratios. It's backed by a very fast entropy stage, provided by Huff0 and FSE library.

We read this blog post on using it to compress fastq files and were intrigued so we ran some tests. We couldn't reproduce the same level of compression that the article claims but that might be due to different input files or the way an index was/wasn't created.
Even in our tests though, we did see that it provides better compression than gzip. Here are the specific values for one of the tests, run using an R1 file of a human WGS run:

Uncompressed fastq
bgzip-compressed
zstd-compressed

179G
58G
47G

This was done using -z 19, to set the compression level to 19 (we found this to be the best tradeoff between size and time to compress). We also saw that using a human genome fasta file to create a compression index, and then using that when compressing human fastq files, actually resulted in worse performance (larger compressed file) than not using the index. But it is entirely possible that we just did it wrong, we didn't spend a lot of time on it.
So you might want to investigate zstd as well, maybe you can get better performance if you play around with the indexing. On the plus side, it is free
and open source. On the other hand, it won't do much with bam or cram files.

Thanks for the input. I'm pretty sure that with Spring, the output size might be around 20-30 GB, if not slightly lesser. Granted, you would have to compress / decompress the data, but the final file size would be pretty small. We have compressed WGS samples with spring, it goes from total 50GB of gz files to around 15-20 GB in final spring file.


text:
compare fasta sequences in pairs and collect metrics
I have 96  fasta files (A1, A2, A3...) from one  plasmid assembly pipeline, and I have another 96 fasta  files (B1, B2, B3 ...) from another plasmid assembly pipeline.
I would like to compare pair wise (A1 to B1, A2 to B2, A3 to B3 etc) and collect some metrics to see how this two assembly method differs.
I am thinking of doing blastn -query A1 -subject B1, but this created a lot visual alignment, rather than detailed align metrics.
Thank you for your suggestions.

OK, I think maybe I will just stay with blastn and do this: blastn -query A1 -subject B1 |grep -E "Score|Identities|Strand" >res_A1_B1
Then I will get report like this  in res_A1_B1, Score = 1.048e+04 bits (5674),  Expect = 0.0
Identities = 5674/5674 (100%), Gaps = 0/5674 (0%)
Strand=Plus/Minus
I can further collect them into tables.

It's easy; just use output format 6, this gives you a nice table that's easy to parse:
blastn -db plasmid2.db -query plasmid1.fasta -evalue 1e-6 -outfmt "6 pident evalue bitscore"

I don't know what parameters you prefer.
... "6 pident evalue bitscore" ...

You just edit the pident, evalue and bitscore above.  You might want to place staxid or something to tell you what it is. A selection of parameters are as follows:

staxid ssciname qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs sstrand

You might want to make the command a bit fancier
blastn -db plasmid2.db -query plasmid1.fasta -evalue 1e-7 -outfmt "6 pident evalue bitscore" | sort | uniq -c | sort -n >> myplasmid.txt

source: Stack Exchange

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

sort: votes


text:
Parsimony-distance neighborhoods in tree space
I have the following question about the paper A parsimony-based metric for phylogenetic trees by V. Moulton and T. Wu.
In this paper, the authors define a "parsimony-distance" between phylogenetic trees on the same set of taxa. They study how this distance is related to other common distances used in phylogenetics, such as counting tree bisection and reconnection (TBR) moves.
For each notion of distance, we can define the "neighborhood" of a tree as the set of trees which have distance one from the first specified tree.
It is claimed in the paper that the parsimony-distance neighborhood contains the TBR-distance neighborhood, for any tree.
Lemma 6.2 states that only caterpillar trees have unit distance to some other tree. This seems to suggest that for any non-caterpillar tree, the parsimony-neighborhood should be empty.
However, the TBR-neighborhood of any tree should be non-empty, since it is possible to apply a TBR move and obtain a different tree.
Are there some assumptions I am missing about when the parsimony-neighborhood is empty, or when the TBR-neighborhood is contained in the parsimony neighborhood?
Also, when applying the formula in Corollary 6.6 to the non-caterpillar tree on 6 leaves, I am unsure how to apply the formula to give the neighborhood size  $|N_p| = 0$, which should be implied by Lemma 6.2.

Another term for "caterpillar tree" is "ladderized tree", the former seems more common in the math community.

The divergence in language makes it difficult to conceptualise and without that it's difficult to approach algebraic abstraction, albeit I've worked on tree space and its landscapes. I suspect I get point 3, but not much more.

"Caterpillar tree" is a term in graph theory. If that is the definition used here then every parsimony tree would be a caterpillar tree, there would never be a "non-caterpillar tree". This is because each branch length defines one "step": "parsimony's version of a mutation", so each "step" within the branch could be considered to be a "branch-less node", or else bifurcating node.

Ladderized tree I don't understand that. Maybe Np = 0 means a complete polytomy, i.e. a star tree, but could be a branch length = 0, i.e. genetically identical taxa.

The concept for an empty neighborhood ... possibly if all topologies are equally parsimonious this could degenerate into an empty tree neighborhood, again a star tree. What parsimony does is always define a bifurcating tree regardless, over replication all alternative solutions may result collectively in an "empty neighborhood". This might be fine from an algebraic perspective, its weird from algorithm perspective though.

When a TBR-neighborhood is equal to a parsimony neighborhood that actually makes sense. I define that as the TBR-neighborhood being "trapped" within the parsimony neighborhood and thats convergence. Here the "step length" has reached its minimum plateau and that tree is the most parsimonious tree ... that makes complete sense. This also makes sense within the context of the paper (I haven't read it - but I can see where they are heading).

What the authors want to do is pick holes in TBR so they want to define the point at which TBR has reached its most parsimonious tree ... I assume they want to show TBR can make errors and their solution is better. Thus they need to define the "end point". At this "end point" they want to show their measure outperforms TBR to reach a lower step size.
Thus they will want to imply that TBR should be replaced with their derivation. Its a nice methodological way to basic (or "fundamental") phylogentics.

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


text:
How much of the space of motifs have been discovered?
I notice when using tools like HOMER there are two different types of outputs: de novo motifs vs known motifs and  it seems no matter what the data is (type of ChIP etc.) there is always a "de novo" motif set in the results.
I was hoping to reach out to experts on insight about how much of the space of motif based structure, within a genome, has been identified. Is there any consensus on what portion of the motif like patterns, for example in human or mouse genome, have been identified? When I find de novo motifs using HOMER in my ChIP-seq or ATAC peaks, what are the chances that in my common-place ChIP experiment (PU.1 for example) I find de novo motifs, not recorded in the databases (i.e. JASPER).

Thanks. I previously closed this question mostly due to the "why" in the title. I've reopened the question because it looks more likely to encourage specific answers.

De novo means that the software infers presence of motifs based on the provided sets of DNA sequences, most commonly peaks from ChIP-seq and similar experiments.
So essentially it stacks the sequences on top of each other and checks for presence of conserved nucleotides that form sequences of certain lengths, for example < 10bp.
Known means that each peak is checked whether it contains a known motif (so is similar to a given known motif) and then in the end a per-motif statistic is calculated to decide whether enriched or not.
If the investigated protein is well-studied then de novo and known is likely similar or identical. De novo does not mean that found motif is novel in terms or unknown or undescribed. It just means that it is derived from the data directly rather than matching known motifs onto the peaks. It's simply two different ways to approach the question. Hope that makes sense.

Thanks @Atpoint, greatly useful and helpful answer.  Part of my question remains, albeit in a HOMER independent context, as to how much of the motif space remains undiscovered, if any?

Probably quite some. Imagine just for how many transcription factors there are no antibodies available.


text:
Singularity & Picard
I am trying to run Picard in Singularity and try as I might I cannot seem to find the path to use. For example with GATK 4.4.0.0 I use the following path:
singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \

where singularity_image is a variable that encodes the path to the gatk.sif. However, I cannot seem to figure out the path I need to use for picard.
When I look inside the Singularity file with the following code:
singularity exec picard_latest.sif /bin/bash -c "ls /"

I get the following:
bin  boot  dev  environment  etc  home  lib  lib32  lib64  libx32  media  mnt  opt  proc  root  run  sbin  singularity  srv  sys  tmp  usr  var

Attempt 1: Using Find
singularity exec picard_latest.sif find / -name "picard.jar" 2>/dev/null

However, this gave no results.
Attempt 2: Going Down The Opt Path
singularity exec picard_latest.sif /bin/bash -c "ls /opt/"

This led to a file called java and that had only one file in it called openjdk so that seemed to be where java was stored.
Attempt 3: Looking In Boot
singularity exec picard_latest.sif /bin/bash -c "ls /boot/"

This yielded no results as it seems to be an empty file.
I have made other attempts but I cannot seem to establish the path I should use. I cannot change the java directly as it is preset on my cluster.
Does anyone know what the path would be for the latest version of Picard or could point me to how to look for this as I cannot seem to find this?

Assuming you've pulled the latest broadinstitute/picard container, your first attempt should have produced the expected result:
$ singularity exec picard_latest.sif find / -name "picard.jar" 2>/dev/null
/usr/picard/picard.jar

And then to run Picard and display the usage, you would use:
$ singularity exec picard_latest.sif java -jar /usr/picard/picard.jar -h
USAGE: PicardCommandLine <program name> [-h]

...

Note that there's now a newer release of Picard, which uses a much newer version of Java. There's no tag for this on Docker Hub (yet). So if you wanted the 'latest', you could try using the quay.io biocontainer:
$ singularity pull docker://quay.io/biocontainers/picard:3.0.0--hdfd78af_1
$ singularity exec picard_3.0.0--hdfd78af_1.sif picard -h
USAGE: PicardCommandLine <program name> [-h]

...

Note that the Bioconda build system also builds Docker containers and these are uploaded to quay.io as part of the Biocontainers project.

singularity exec picard_latest.sif /bin/bash -c "ls /usr/picard/" does show the jar file there however, singularity exec picard_latest.sif find / -name "picard.jar" 2>/dev/null does not seem to work. Out of curiosity when I use your usr/picard/picard.jar path I get an out file with an error Error: Invalid or corrupt jarfile /path/to/picard_latest.sif - the path is correct why would it give this error? The main two lines of my code are java -jar ${dir} /usr/picard/picard.jar \     BedToIntervalList \

@Indira Do you mean where dir is a path to your picard_latest.sif? I think you may have pointed to your GATK sif file. When you say the find command does not seem to work, what do you mean exactly? If it is just slow, try find /usr/picard ...

I wrote it wrong in the comment - the path for the dir is to the picard_latest.sif

When I say the find command did not work I mean that all it did was not produce a result - it kept running but didn't produce any result.

@Indira So did the command complete with a zero exit status or did you end up needing to interrupt it with ? I see the mistake in your code, it should be: singularity exec "${sif}" java -jar /usr/picard/picard.jar BedToIntervalList

Thank you. The above worked. It seems like I have multiple small errors but they should be fixable.


text:
Integrating bulk RNA-Seq data with different sequencing depths and from different sources
I am attempting to integrate different bulk RNA-Seq datasets. While this is not ideal, I'm trying to reduce the technical variability in these datasets by using data generated by similar protocols (for example, all of them have been sequenced with a Poly-A enrichment strategy and Illumina short read sequencing). I am aware the results I obtain from this analysis will be very limited, but still want to pursue the analysis. I'm looking for advice on two aspects of the analysis, but first some more info about the data:

There are 5 bulk RNA-Seq datasets and 4 biological groups. Each datasets only comprises one biological group. Only two biological groups have two datasets, and such datasets overlap in PCA.
In Principal Component Analysis (PCA) and the first PC (~26% of variation) is very correlated with the biological variable (cell differentiation), suggesting there is a biological signal among all the noise/technical bias.

Despite these preliminary results I'm not sure I'm accounting correctly for some aspects: The datasets have been sequenced with different sequencing coverage (in fact, different numbers and length of reads). Some are single read and some paired-end. My questions are:

For PCA, I'm normalizing with variance stabilizing transformation (VST). How should I account for the differences in read numbers and type between datasets?

For differential expression (DE) analysis I'm using DESeq2. How should I account for the differences in read numbers and type between datasets? There is no batch variable, as all datasets come from different sources (batches, laboratories)

I'm happy to provide more information if necessary. Thanks

Welcome to the site. I'll provide an alternative view to @gringer's early next week, the results may be identical however. I agree with the theme of the post and the encouraging results from PCA. What I'm going to suggest is leveraging PCA as a diagnostic on the standardisation for downstream analysis

Per-sample coverage / count differences are expected in all differential expression analysis, and correction for that is built into DESeq2. You shouldn't need to do any additional correction with respect to proportional counts when looking at a PCA derived from the VST.
For accounting for different sources / types / etc, it's a good idea to add that information into the metadata, and test for differential expression with and without accounting for those variables to see if anything changes. This is recommended in the DESeq2 manual:

If there is unwanted variation present in the data (e.g. batch effects) it is always recommend to correct for this, which can be accommodated in DESeq2 by including in the design any known batch variables or by using functions/packages such as svaseq in sva (Leek 2014) or the RUV functions in RUVSeq (Risso et al. 2014) to estimate variables that capture the unwanted variation. In addition, the ashr developers have a specific method for accounting for unwanted variation in combination with ashr (Gerard and Stephens 2017).

See the multi-factor design section of the manual for more information (which has a "single-read" vs "paired-end" example):
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs
If after accounting for the different variables there is not enough information to properly separate groups, DESeq should tell you. Ideally you should aim for six replicates of each state of each independent variable (easier to do when there are interacting variables), but it's quite rare to see that in differential expression analysis.

"There are 5 bulk RNA-Seq datasets and 4 biological groups. Each datasets only comprises one biological group. "-- to me it seems the data is confounded and it may be impossible to disentangle biological and technical differences

In case of different sequencing depth one performs rarefaction, to make the data comparable.
The standard approach is rarefaction by size, when one takes the number of reads in the smallest sample, and samples without replacement the same number of reads from the other samples. Potentially, one could carry out upsizing to a greater number of reads (e.g., if the smallest sample is an outlier in terms of the number of reads, and downsizing may result in losing too much information.) I suggest reading the tutorial for MOMR package for more detailed description (the R code is also available.)
There are more sophisticated methods of rarefaction, such as rarefaction by coverage - these are of significant scientific interests, but tricky to use and not very widespread. See, e.g., this article (it deals with sampling species, but generalization to reads is straightforward... though reads are not species.)
Remark: Of course many good packages (like those mentioned by @gringer) might already perform rarefaction by default - but it is always a good idea to check the documentation.


text:
No common application protocol between the client and server? An error when using Entrez Direct
After a long time not running EDirect, I got an error when I tried using it. It said," curl: schannel: next InitializationSecurityContext failed: Unknown error (0×80090367) - There is no common application protocol between the client and server...
I tried using WiFi and hotspot, but it did not work.
Additional info:
1- I ran EDirect on Cygwin64 Terminal on my personal computer with Windows 11 operating system.
2- I first ran 'export PATH=${PATH}:${HOME}/edirect'. After that, I tried 'esearch -db pubmed -query "stroke AND english [FILT] AND has abstract [FILT] NOT review [PTYP]" | efilter -mindate 2000 | efetch -format xml | xtract -pattern PubmedArticle -element PMID AbstractText > dm_cs.txt'
I used to ran it successfully on the same PC over WiFi or hotspot three month ago.

Please [edit] your question and add: i) your operating system, ii) the exact command you tried to run and iii) tell us if there is anything special about your setup (are you behind a firewall? are you connecting over ssh or something similar? anything else relevant).

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

One possibility is it's simply NCBI blocking the request to avoid being overrun by e.g. bots. This occurs if you've been running the script repeatedly within a single day. It will be fine tomorrow in this scenario.
The second possibility is you've not supplied an email address.
The final possibility is if you're running a VPN. NCBI has taken a strong dislike to VPN's per se, I've never figured out its algorithms - it seems to permit a fixed number of requests then pointedly refuses. Just turn off the VPN.

If its point 2 then understanding how you've coded this would be useful to allow an email address to be sent. Biopython has a function for this. I just use Biopython BTW.


text:
calculate mismatch frequency/rate from a BAM file
I am working with PAR-CLIP data where experimental procedure induces T to C mismatches. The other types of mismatches we can consider coming from artifact such as PCR errors. It could also be an existing SNP. For simplicity, let's consider T to C mismatches are experimentally induced and rest of the mismatches are background mismatch rate.
From Illumina NGS seq after performing alignment, for each of he read groups (reads collapsed in a 30-70 base region with minimum 10 read coverage). Up to this point its already done.
I want to compare T to C mismatches vs background mismatches. To do so I need to calculate the rate of signal mismatches (T to C) and background mismatches.
Using samtools mpileup I can get total number of all kind of mismatches for a read group, however, takes some time.

How do I calculate background mismatch rate from this?
How do I calculate total mismatches in a bam?

I was thinking about something like
for a read group:
mismatch_rate_group=total_mismatches_group/(total_read_group*total_base_group)
for a BAM file:
mismatch_rate_global=total_mismatches_bam/(total_reads_aligned_bam*genome_size_base)

What important factor I might missing in my approach?

I  am not using any variant caller (e.g. bcftools). Initially I tried with samtools mpileup and then used a shell script to figure out the mismatches. But I am not sure is this is the right thing to do. Or even there are any easier methods. If one read has a mutation, I consider that a mismatch and I presume it to be artifact and/or PCR errors.

Please [edit] your question and explain this. You really should use a dedicated variant caller. NGS is error prone, and you cannot just assume that >1 read means the variant is real. In any case, writing your own thing to do this would be far more work and would almost certainly result in worse data than if you use a dedicated variant caller which has been designed for the job. If [edit] your question and give some more details (what sequencing technology, what species, what coverage etc) we will be able to recommend some tools.

@terdon , thanks for your suggestion. I tried to add some more details in my question. I hope this is better now.
Additionally, NGS is error prone : I am less interested in the variant, rather I want to know the rate of such errors so that I can differentiate my signal from all of such errors.

I am not using any variant caller

Use a variant caller. bcftools, at least, will output variant frequency per variant:
bcftools mpileup -ugf ref.fa sample.bam | bcftools call -mv > output.vcf

As an alternative approach, I wrote some code to parse the mpileup output and extract counts for different alleles; this output may be easier for you to parse and calculate your desired frequencies:
# readstomper.pl -- generate base-packed proportion statistics for
#   output from 'samtools mpileup'
#
# example use:
# samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
#   /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.txt.gz
#
# example output:
# Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
# Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
# Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
# Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
# Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
# Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
# Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
# Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00


text:
Manta Evidence Folder Is Empty
I am looking at Manta on Github and I ran a bit of code:
The following variables were defined in the code:

input_dir is the file path to a folder with less than 10 .bams as per the limitations of Manta
output_dir is the path to the output
MANTA_INSTALL_PATH is the path to the Manta installation
reference is the path to the reference fasta

    for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
      mkdir -p ${sub_dir}

      # Use reference fasta path specified by user
      ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
        --referenceFasta ${reference} \
    --runDir ${sub_dir} \
        --exome
    done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

The run worked and each sample had a folder with the following files:
results  runWorkflow.py  runWorkflow.py.config.pickle  workflow.error.log.txt  workflow.exitcode.txt  workflow.warning.log.txt  workspace

Each workflow.error.log.txt had nothing in it implying the the run had been successful. Each results folder had 3 subfolders
evidence  stats  variants

The stats and variants folder each had outputs but the evidence folder had nothing. Is this normal? Should it have something in this folder? Why is this folder there in Manta if nothing is in it?

The documentation on the GitHub page you linked to has this to say about the evidence folder:

Generating evidence bams
Using the --generateEvidenceBam option, Manta can be configured to
generate bam files of evidence reads for SVs listed in the candidate
vcf file.
It is recommended to use this option together with the --region
option, so that the analysis is limited to relatively small genomic
regions for debugging purposes.
The evidence bam files are provided in
${MANTA_ANALYSIS_PATH}/results/evidence, with a naming format
evidence_*.*.bam. There is one such file for each input bam of the
analysis, containing evidence reads of the candidate SVs identified
from that input bam. Each read in an evidence bam keeps all
information from the original bam, and it contains also a customized
tag in the format:
ZM:Z:${MANTA_SV_ID_1}|${EVIDENCE_TYPE},${MANTA_SV_ID_2}|${EVIDENCE_TYPE}.
For example, ZM:Z:MantaINV:5:0:1:0:0:0|PR|SRM,MantaDEL:5:1:2:0:0:0|SR
One read can have more than one of the three evidence types: PR for
paired reads, SR for split reads, and SRM for split read mates. One
read can be evidence for multiple SVs, which are separated by commas
in the tag. Notice that the number of evidence reads for a particular
SV in the evidence bam files could be more than the evidence counts
(PR and SR) in the final vcf files. This is because more stringent
criteria are applied for generating evidence counts in the final vcf
files.

So it looks to me like you need the --generateEvidenceBam option to generate evidence, but that the directory evidence is created regardless of whether you use the option. Try running the analysis again, but using --generateEvidenceBam this time.

Thank you I checked the documentation again this seems to be the case. It is just an empty folder until --generateEvidenceBam is specified.


text:
Strategy for merging many VCFs
Can anyone suggest a strategy for speeding up VCF merging?
I have ~44,000 single-sample VCFs that I am trying to merge into a multi-sample VCF with bcftools merge, but the job keeps timing out on the cluster (just timed out after 4 days). I could obviously try running it for a couple of weeks, but I was hoping there is a better strategy someone can suggest?
The exact command I am running is
bcftools merge --no-version -m all -i 'DP:avg,RO:avg,AO:avg' -F + -o merged.vcf --file-list merge.fofn

Happy to try any other tool that will give the same merge results.

We created this over the pandemic when we needed to merge loads of VCFs which all have the same header and first few columns (ie same variants)
https://github.com/iqbal-lab-org/ivcfmerge
That page gives detailed instructions but i'll paste a couple of examples
of how you can use it
== In Python
=4.1.1 If the number of input files is small (can be opened all at once)
from contextlib import ExitStack
from ivcfmerge import ivcfmerge

filenames = [...]    # List/iterator of relative/absolute paths to input files
output_path = '...'  # Where to write the merged VCF to

with ExitStack() as stack:
    files = map(lambda fname: stack.enter_context(open(fname)), filenames)
    with open(output_path) as outfile:
        ivcfmerge(files, outfile)

=4.1.2 If the number of input files is big (cannot be opened all at once)
from ivcfmerge import ivcfmerge_batch
filenames = [...]    # List/iterator of relative/absolute paths to input files
output_path = '...'  # Where to write the merged VCF to
batch_size = 1000    # How many files to open and merge at once

ivcfmerge_batch(filenames, output_path, batch_size)

Pretty sure we initially did batched bcftools merge, I can't honestly remember why we moved to this. Thankfully not having to do these merges very often

Thanks Zam. Unfortunately I fail the first two assumptions listed in the README.

I ended up using an iterative approach (inspired by this post), which took 1.5 hours - so a massive speed up (and that was just the sequential version).
Note: fd is analogous to find, but way nicer to use.
First, merge subsets of n VCFs
# size of the subsets
n=200
# create file of filenames for each subset
fd -e vcf . path/to/vcfs/ | split -l $n --additional-suffix '.fofn' - subset_vcfs

-e vcf tells fd to find all files with a .vcf extension, you can obviously change this to .bcf/.vcf.gz etc. if necessary. split splits the list of all VCFs files into subsets of size -l (200) and creates the file of filenames with the prefix subset_vcfs. It will add a random suffix to this - e.g. subset_vcfsav. Plus the additional suffix .fofn. See the split manual for all options.
Now iterate through these subsets and merge the VCFs listed in them.
Subset merge option 1 - sequential
This just merges the VCFs one at a time
for sub in *.fofn;
do
    echo "Merging subset ${sub}..."
    out_vcf="merge.${sub}.bcf"
    bcftools merge --file-list "$sub" -o "$out_vcf"
    bcftools index "$out_vcf"
done

Subset merge option 2 - parallel
You can easily parallelise the previous step using fd on the subset files. See this post for more details. Or basic bash for loop parallelisation.
Create a script (merge_subset.sh) of what we want to run on each subset
#!/usr/bin/env bash
sub="$1"
echo "Merging subset ${sub}..."
out_vcf="merge.${sub}.bcf"
bcftools merge --file-list "$sub" -o "$out_vcf"
bcftools index "$out_vcf"

Now we run this script with each subset file of filenames using fd
# number of processes to use. remove this option to use all available
jobs=4
fd -j $jobs -e fofn -x bash merge_subset.sh '{}'

Final merge
Now we merge all of these VCFs into our final VCF
ls merge.*.bcf > merged.fofn
bcftools merge -o merged.bcf --file-list merged.fofn

Perform a tree based merge of the files.

Sort the files into groups that can easily be merged in a measurable amount of time.
Merge each group into a separate Variant Call Format (VCF) file

Has the benefit that files can also be grouped by type / source / run
Further benefit that small merges can be checked to ensure you like the output before proceeding to the next stage.
Further benefit that there are also subgroup files available to be distributed or shared separately.
Downside, extra storage space until you have a final product.

Merge subfiles upward in stages until you arrive at a complete 44,000 entry file

If you're happy with the earlier test results, process can be automated with relatively simple iterative shell scripts.
Possibly faster, as the bulk merge needs to load every single file into memory, and you may be hitting a capacity bottleneck.
Worked in Fluid Dynamics for years, and often ran into a similar issue.  Could run many small cases representing a large case faster than running the large case, because large case would overwhelm memory and transfer capacity.
In some clusters you may also be causing nodes to silently crash, which then stalls the entire process without a visible notification.


text:
Custom-made kraken2 catalogue
I am looking for instructions on creating a custom-made kraken2 catalogue from the sequences of my choice and the corresponding phylogenetic tree in newick notation.
I will appreciate an input on whether/how it can be done. The prescription provided in the kraken2 manual seems to assume that sequences and phylogeny originate from the NCBI. Struo2 extends this to GTDB sequences. However, in my case the sequences are collections of concatenated genes (i.e., not full genomes), provided by a third party.
Remark: (in response to comments) I don't need to work with newick format - I have the tree in newick format. So converting it to dmp might be a part of the answer.

Newick has a lot of different ways to represent data, so it's helpful to know specifically what your input data looks like. When asking a question, it's helpful to provide both example input and example output. You've provided good information on the output side through the link, but no information on the input format.

@gringer Does kraken2 need this information? The question is not about solving my particular research problem (which is off topic?), but about creating a kraken2 catalogue from sequences and phylogeny, other than those readily provided by NCBI or GTDB (via Struo). [contd.]

You seem to be suggesting, that, once phylogeny is converted to dmp format, the problem is easily resolved - if that's the case, then state this explicitly, specify the tools for converting newick to dmp, and the minimal information that a newick tree should contain for creating a meaningful kraken2 catalogue. The would make for a nice answer, and useful for other users that stumble on the same problem.

I'm not trying to give answers, just clarify the question. My comments are more generally applicable to all questions on this site: a good question for StackExchange has a specific, well-defined input, and a well-defined output. All I'm trying to do is make sure that this question will receive good answers.

Yes, phylogeny is important. Kraken2 is a phylogenetic mapper, and uses the tree structure to work out the best match when it cannot be assigned reliably to a single species.


text:
Can bwameth accept unaligned Bam?
Does bwameth accept bam (unaligned) type as input or I will have to convert unaligned bam to fastq (samtools fastq bam) to use it as input?

No, it's a bwa mem wrapper for bisulfite data essentially and bwa mem accepts fastq only.

You can pipe uBAM to bwamem with samtools fastq in.bam | bwa mem -p - > out.sam, provided that paired-end reads are grouped by name.


text:
What does the IGV refseq genes with three different thicknesses of lines mean?
What does the IGV refseq genes with three different thicknesses of lines mean?

In gene tracks, exons ( regions of the transcript that aligning with the genome ) would be depicted as blocks, while introns would be drawn as the horizontal lines that connect the exons.  The direction of transcription would be indicated by arrowheads on the introns.  Coding regions of exons would be depicted as tall blocks, while non-coding exons are shorter.  The 5'-most exon of each transcript is shorter on the left, indicating an untranslated region, and is taller on the right, indicating a coding sequence.

The thinkness indicates whether something is an intron, a coding exon, or a UTR. The very thin lines are introns, the full thickness boxes are coding exons and the thinner boxes are UTRs.
UTRs, untranslated regions, are exons, they are included in the final mRNA, but they don't code for protein. So, by convention, they are shown as thinner exons to give the user a visual indication that although this is still an exonic region, it's a special one.
I don't really know where this is documented, it's pretty standard across all genome browsers and bam visualizers I know though so it's a common idiom.


text:
DESeq2 throwing error while normalizing raw microarray expression data due to presence of negative values
This question was also asked on Biostars
I am trying to download and analyze a miRNA expression dataset from NCBI GEO (GSE25631). I specifically want non-normalized data to perform normalization and my own set of other analyses later on. Accordingly, I downloaded the Series Matrix File (to match the Sample IDs with their phenotype) and the non-normalized data GSE25631_non-normalized_data.txt.gz programmatically using the GEOquery package in R. To normalize the data I'm using the DESeq2 package. When I run the DESeqDataSetFromMatrix function I get this error:

Error in DESeqDataSet(se, design = design, ignoreRank) :
some values in assay are negative

I went back to the non-normalized data and observed that there are a lot of negative values. I need suggestions on how to deal with this issue. Should I replace all such values with NA? I also observed that every probe has a detection p-value associated with it. Should I also filter out the expression values which have a p-value higher than 0.05?
> head(raw_counts)
        gbm-0001 Detection Pval gbm-0011 Detection Pval gbm-0013 Detection Pval gbm-0024 Detection Pval gbm-0058 Detection Pval
        ILMN_3167304   27.875      0.3071961 -24.9375      0.6689599    14.25       0.350027   4.3125      0.4723163     8.25      0.3966338
        ILMN_3168038  211.875       6.42E-05   1.0625      0.4925717    99.25    0.003645935  60.3125      0.1657064    72.25     0.01086376
        ILMN_3167890 13839.88       3.68E-38 23507.06       3.68E-38 22345.25       3.68E-38 20482.31       3.68E-38 21419.25       3.68E-38
        ILMN_3167526    5.875      0.4577177 105.0625     0.03279042     8.25      0.4117529 -29.6875       0.683706    11.25      0.3604082
        ILMN_3167209   82.875     0.06708142 182.0625      0.0007096    60.25     0.05167192   9.3125      0.4403947    72.25     0.01086376

Here's the R script which I used for all the mentioned steps.
library(GEOquery)
library(DESeq2)
library(tidyverse)

data <- getGEO(GEO = "GSE25631",
           destdir = "E:\\GSE25631",
           GSEMatrix = TRUE,
           AnnotGPL = FALSE,
           getGPL = FALSE,
           parseCharacteristics = TRUE)

clindata <- data$GSE25631_series_matrix.txt.gz@phenoData@data

url = "https://www.ncbi.nlm.nih.gov/geo/download/?
acc=GSE25631&format=file&file=GSE25631%5Fnon%2Dnormalized%5Fdata%2Etxt%2Egz"
download.file(url, "E:\\GSE25631\\GSE25631_raw_counts.gz")
raw_counts <-read.delim("E:\\GSE25631\\GSE25631_raw_counts.gz",
                     stringsAsFactors = FALSE, sep = "\t")
colnames(raw_counts) <- raw_counts[5,]
raw_counts <- raw_counts[-(1:5),]
raw_counts <- raw_counts %>% remove_rownames %>% column_to_rownames(var = "ID_REF")
raw_counts <- raw_counts[ , c(TRUE, FALSE)]

colnames(raw_counts) <- clindata$title
rownames(clindata) <- clindata$title
colnames(clindata)[colnames(clindata) == "disease state:ch1"] <- "Condition"
clindata$Condition[clindata$Condition == "primary glioblastoma multiform (brain tumor)"] <- "GBM"
clindata$Condition[clindata$Condition == "normal"] <- "Normal"
clindata$Condition <- as.factor(clindata$Condition)

dds <- DESeqDataSetFromMatrix(countData = raw_counts,
                          colData = clindata,
                          design = ~ Condition)

Answered at biostars: https://www.biostars.org/p/9561238/#9561240

DESeq2 is meant for digital counts (RNA-seq and similar assays). Arrays are not supported, nor is the implemented method meaningful for array intensity values.
Guides on how to handle these data are probably in here:

Thanks for linking to the other answer, but when you do so, please also reproduce the answer here so that it will not depend on the external link not breaking.


text:
How do i search by these gene properties in NCBI's GenBank?
I want to search GenBank entries by gene product and taxon id.
Taking the 16s rRNA as an example for a gene and Piscirickettsia salmonis (taxon:1238) for a taxon:

https://www.ncbi.nlm.nih.gov/nuccore/219846389

How do I retrieve all GenBank entries that have the following properties /product="16S ribosomal RNA" and /db_xref="taxon:1238"?
Update: It would be appreciated if you also showed me how to retrieve GenBank entries like this, which don't have an exact match to the "16S ribosomal RNA" product.

The research study that submitted this sequence was in chile, but I think I've seen other sequences for the same organism, also obtained from chile.

"The salmon, which displayed typical disease signs, was collected from a seawater netpen in southern Chile where an epizootic was in progress. "
https://www.microbiologyresearch.org/content/journal/ijsem/10.1099/00207713-42-1-120

This can be done using the NCBI eutils command-line tools. For taxa, the correct incantation for the search string is txid<taxID>[Organism] (assuming you want that taxon and all its children), but I haven't found the right one for the product field. However, it can be approximated by searching for a string in the tile line via <search>[TITL]. This seems to work fairly well, assuming you want the fasta files (the -format argument can be changed to genbank if you want the genbank record instead):
$ esearch -db nuccore -query "txid1238[Organism] AND 16s[TITL]" | efetch -format fasta | grep -A 1 '^>' | head
>MW433822.1 Piscirickettsia salmonis strain HRV16 16S ribosomal RNA gene, partial sequence
TCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAACATGACTAGCTTGCTAGTTGA
--
>MW433821.1 Piscirickettsia salmonis strain HRV13 16S ribosomal RNA gene, partial sequence
AATTTCAGGTGTAGCGGTGAAATGCGTAGATATCTGAAGGAACACCGGTGGCGAAGGCGACTTTCTGGAT
--
>MN023100.1 Piscirickettsia salmonis clone Ps18621-6 16S ribosomal RNA gene, partial sequence
GCTTGCTAGATGCTGACGAGTGGCGGACGGGTGAGTAACGCGTGGGAATTTACCTTTTAGTGGGGGATAA
--
>MN023099.1 Piscirickettsia salmonis clone Ps18429-3 16S ribosomal RNA gene, partial sequence

Unfortunately I can't find anything for searching specifically in the features of an NCBI record. There'll be some way to do it; I just don't know how.
Based on the updated request, it's possible to do a full-text search through the entire metadata (e.g. for 16S ribosomal RNA) by adding the search term without specifying a field. Depending on what you want to do with the results, and how many were produced, it could be easy enough to manually unwanted sequences (e.g. whole genome shotgun assemblies, or chimeric products):
$ esearch -db nuccore -query "txid1238[Organism:exp] AND (16S ribosomal RNA)" | efetch -format fasta | grep -A 1 '^>' --no-group-separator | tail -n 12
>EU170254.1 Piscirickettsia salmonis strain NVI 5692 16S ribosomal RNA gene, partial sequence
CACATGCAAGTCGAACGGTAGCAGATCTAGCTTGCTAGATGCTGACGAGTGGCGGACGGGTGAGTAACGC
>U36941.1 Piscirickettsia salmonis 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCCTGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U55015.1 Piscirickettsia salmonis 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCCTGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U36942.1 Piscirickettsia salmonis strain NOR-92 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCCTGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U36940.1 Piscirickettsia salmonis strain EM-90 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCATGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U36915.1 Piscirickettsia salmonis 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCATGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT

FWIW, here's the list of all the searchable fields in the nuccore database:
$ einfo -db nuccore -fields
ACCN    Accession
ALL     All Fields
ASSM    Assembly
AUTH    Author
BIOS    BioSample
BRD     Breed
COMP    Component Accession
CULT    Cultivar
DIV     Division
ECNO    EC/RN Number
FILT    Filter
FKEY    Feature key
GENE    Gene Name
GPRJ    BioProject
ISOL    Isolate
ISS     Issue
JOUR    Journal
KYWD    Keyword
MDAT    Modification Date
ORGN    Organism
PACC    Primary Accession
PAGE    Page Number
PDAT    Publication Date
PORG    Primary Organism
PROP    Properties
PROT    Protein Name
SLEN    Sequence Length
SQID    SeqID String
STRN    Strain
SUBS    Substance Name
TITL    Title
UID     UID
VOL     Volume
WORD    Text Word

To specifically answer the question the search that is used is
((1000:1500[Sequence Length]) AND Piscirickettsia salmonis[Organism]) AND 16S

That is place here and select 'nucleotide'.
Genbank uses taxon id internally, it's a database reference not a biological code, and can't be used as input - its completely unique to Genbank. Anyway the search results in 276 records, and are shown here, this is near full length 16S thats cool.
Clarification [taxon id] is not a field in Genbank, it's a database code, the tag is [organism].  Genbank knows what 16S is - its the most famous genetic marker in microbiology - so there is no need for a [tag]
At the base of the page there is a "Send to:" option (right hand corner). Click and then select "File". The 276 sequences will then download.
If [Sequence length] is removed this result 3820 sequences here. There are some 16S-23S sequences on the database result in a 2,198 bp length, so extending the sequence length will trap these by 1000:2200[Sequence Length]. There are 431 bp sequences - which are not cool at all, probably a lot.
The approach used below identifies nearly 1000 is complicated it traps everything of relevance which is close to full length and cuts out the 16S from  2,198 bp length and beyond. The 431 bp partial 16S sequences have been removed. Thus the numbers below are the idealised data set, so its about 1000 sequences and the tree will guide you into major groups - which is probably what is wanted, i.e. to place a 16S from a local salmon farm into the broader clade diversity of the species.

The answer is for this species is there are around a thousand sequences on Genbank for the species of interest which are split into 6 identical genotypic groups. The demonstration of the species is represented by the scale bar, i.e. it is very little genetic diversity.

Cycloclausticus sp is the outgroup and each member of the in-group is Piscirickettsia salmonis

Just to address the comments

"Cycloclasticus"? The question is not specific to /db_xref="taxon:1238"

It is specific there's a fundamental principle of trees. A tree needs a root, all the ingroups are all Piscirickettsia salmonis. Without a root, i.e. something that is not Piscirickettsia salmonis there is no direction in the tree. Cycloclasticus is a good or reasonable root.

how do i access this interface?

Its complicated to explain and its a different question.


text:
GDC API: Obtain the Objective Power (or Magnification) of a WSI file
I've been trying to find out a way to retrieve the Objective Power (or Magnification) of a Whole Slide Image file without downloading it first.
I've seen that some WSI's Objective power might be 40x and others 20x, and I need a way to filter the ones I want to retrieve according with this parameter.
However, I cannot figure out the query I need for this.
Can anyone help me
Some details:
I'm using the GDC API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api) to retrieve Whole Slide Image files (in this case .svs files).
I want to retrieve this files according with its Objective Power, since images in the TCGA dataset (https://www.genome.gov/Funded-Programs-Projects/Cancer-Genome-Atlas) might have 20x or 40x of objective power.
One way to do this is to retrieve the whole image using the data endpoint of the GDC API and then use a library like Openslide to check the objective power:
import requests
from io import BytesIO
import openslide

# Set up the GDC API endpoint
endpoint = 'https://api.gdc.cancer.gov/data'

# Set up the query parameters
query = {
    'ids': ['my_file_id'],
}

# Send the request to the GDC API
response = requests.get(endpoint, params=query)

# Retrieve the binary data from the response
binary_data = BytesIO(response.content)

# Open the WSI file using the Openslide library
slide = openslide.open_slide(binary_data)

# Retrieve the Objective Power of the WSI
objective_power = slide.properties.get(openslide.PROPERTY_NAME_OBJECTIVE_POWER)

print('Objective Power:', objective_power)

However, retrieving the whole file takes too much time, so I would like to know if anybody knows of a field in the GDC API that allows me to check the objective power of a WSI without retrieving the file itself.


text:
How do I save selected sequences in seqkit to a file?
Consider:

seqkit grep -nrp ";s__nosocomialis;"  ribogrove_11.217_sequences.fasta

Using this command lets me select all sequences related to the organism "A. nosocomialis" in the RiboGrove database.
How do I save these selected sequences to a new ".fasta" file?

Based on the docs, seqkit grep has a -o "output file" option, so if you get the correct output with:
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta

Then this command will (in theory) give you your desired outcome:
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta -o new.fasta

Thanks, it worked. I don't know how i didn't see the output command while surfing the documentation.

For any program that prints output to the terminal (standard output), you can always redirect to a file using > or >>. Either of the two will create the output file if it doesn't exist, but > will overwrite any contents of an already existing file, while >> will append to it without overwriting. So you could just do:
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta > out.fa


text:
Merging multiple seurat datalist together
I have a R script that subsets out the tumor cells from 6 different individual seurat objects.
It gets a list of file paths for all RDS files containing the expression data for multiple datasets in a directory.
It creates two empty lists, datasets and dataset_index, to hold the subsets of tumor cells from each dataset and their corresponding indices, respectively.
It loops through each file path, reads the RDS file, and subsets the tumor cells from the expression data.
If the dataset does not contain any tumor cells, it prints a message indicating that fact.
The subsets of tumor cells are appended to the datasets list and their corresponding indices are appended to the dataset_index list.
scfiles = list.files('/LungDataSets/',pattern = '.*SCTransform_Symphony.*.RDS',recursive=TRUE)

scfiles = paste('/LungDataSets/',scfiles,sep='')

datasets = list()
dataset_index = list()

for (i in 1:length(scfiles)) {
    scData = readRDS(scfiles[i])

    if (is.element('Tumor', unique(scData$cell_ann)))  {
        tumorData = subset(scData,subset = cell_ann=='Tumor')
        datasets = append(datasets,tumorData)
        dataset_index = append(dataset_index,toString(i))
    }
    else {
        print(c(scfiles[i], 'does not contain Tumor cells'))
    }

    rm(scData)

}

Now I want to merge all those tumor cells from the 6 datasets into one seurat object.
I've tried merge and rbind but they don't work.
megaObject <- merge(datasets)
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("Seurat", package = "SeuratObject")’ to a data.frame

The first parameter of merge should be a Seurat object, the second (y) can be one Seurat object or a list of several. Try:
merge(x = datasets[[1]], y = datasets[-1])

See the merge vignette for more details. You may want to use the add.cell.ids option to be able to tell which dataset each cell originated from.

You should be able to use something like this:
Reduce(
    f = function(x, y) {merge(x, y, merge.data = FALSE)},
    x = datasets # list of Seurat objects
)

This will create a new Seurat object based on the multiple seurat objects in your list. It might be good idea to store the "sample" information within the metadata slots of individual objects. Probably a better approach is to have cell barcodes of each sample prefixed with the sample id.

I can't reproduce things quickly without a minimum reproducible example (e.g. with sample data from the package). However, the error is telling you that R cannot convert the object to a dataframe when trying to use merge() (at least I suspect that this error occurred when you ran merge() since you posted it right under it). I suggest you convert everything you need to a dataframe (there are quite some questions, e.g., on biostars regarding this conversion). Everything you need because a SeuratObject cannot simply be cast to a dataframe, only parts of it, e.g., data_transcripts_seurat@assays$SCT@counts %>% as.matrix() %>% t() %>% as.data.frame() (from this biostars question). However, you might need some other information. Therefore, I suggest you explore the SeuratObject using RStudio and check what information you want to draw from it. After the (interesting) parts of the SeuratObject have been cast to dataframes, e.g., in the for-loop, you can go ahead with merging.

There is no need to convert the Seurat object into a data frame, because a merge() function exists for Seurat objects.


text:
illuminaHumanv2.db returning NA for all Illumina probe ID to gene symbol conversion
This question has been cross-posted on Biostars but has not received any answers yet, so I'm posting it here. I'll update here if it gets answered there.
I have downloaded a miRNA expression dataset from NCBI GEO (GSE25631) to study differential gene expression and perform other analyses. As mentioned in GEO, this profiling was performed on GPL8179 Illumina Human v2 MicroRNA expression beadchip. Accordingly, I installed the illuminaHumanv2.db package to annotate the Illumina probe IDs. However, when I tried to convert the Illumina probe IDs to their corresponding gene symbols, it returned NA for all the probes. Can anybody identify the issue with this approach? Here's the R script which I used.
library(GEOquery)
library(limma)
library(illuminaHumanv2.db)

Sys.setenv(VROOM_CONNECTION_SIZE = 256000)

data <- getGEO(GEO = "GSE25631",
           destdir = "E:\\GSE25631",
           GSEMatrix = TRUE,
           AnnotGPL = FALSE,
           getGPL = FALSE,
           parseCharacteristics = TRUE)

data <- data$GSE25631_series_matrix.txt.gz

raw_intensity <- exprs(data)

samples <- as.character(pData(data)[,"title"])

raw_intensity_log <- log2(raw_intensity)
colnames(raw_intensity_log) <- c(rep("GBM", 82), rep("Normal", 5))

probe_id <- rownames(raw_intensity_log)
mapping <- data.frame(Gene = unlist(mget(x = probe_id, envir = illuminaHumanv2SYMBOL, ifnotfound = NA)))

head(mapping, 10)
               Gene
ILMN_3166935   NA
ILMN_3166938   NA
ILMN_3166940   NA
ILMN_3166941   NA
ILMN_3166943   NA
ILMN_3166944   NA
ILMN_3166945   NA
ILMN_3166948   NA
ILMN_3166952   NA
ILMN_3166953   NA

The GSE25631_RAW.tar file in the supplementary section has all the Probe ID to Gene mappings but I'm curious to know why this approach didn't work.

Have you confirmed that the probe names are valid using a different tool? Could this be caused by having windows line endings in your file and those are included when querying? Or, since your E:\\ path suggests you are using Windows, could it be the other way around and something is inserting a \r before every trailing \n?

I checked the probe-to-gene mappings present in illuminav2SYMBOL. Surprisingly, none of my probe IDs are present, so it returned NAs. The list doesn't seem to contain any mappings for miRNAs. I ended up using the mapping file provided by the manufacturer.

I closed this as non-reproducible since the issue seems to be very specific to your case and your probe names.


text:
Do you need galaxy.yml when you start creating galaxy?
First question: I found that there is no galaxy.yml in the config directory of the running version of docker galaxy, only galaxy.yml.sample file and galaxy.yml.interactivatetools file, so do I need to create it, or do I need to create a new one using the cat command? galaxy.yml file?
The second question is about these initial notes, where can I go deeper?

What have you already tried? Does Galaxy work for you if you run it without a .yml file? Have you looked at the galaxy documentation on the Galaxy website?


text:
GATK SplitIntervals: A USER ERROR has occurred: but no positional argument is defined for this tool
I am trying to create a GATK g CNV pipeline and I need to make shards before running the GermlineCNVCaller therefore, I believe I need to use SplitIntervals.
I have the following predefined variables:

singularity_image is the path to the GATK sif file
bam_file is a bam file from a defined bam directory
reference is a path to a .fa file
targets_interval_list is a interval list created by Picard (ending with .interval_list)
NUM_BINS is a manually preset value from 100 to 200.

  singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \
  SplitIntervals \
  -I ${bam_file} \
  -R ${reference} \
  -L ${targets_interval_list} \
  --scatter-count ${NUM_BINS} \
  -O ${output_dir}/Interval_list_${bam_filename}_.shard_\%${NUM_BINS}\%_\%{shard_index}\%.interval_list

I have tried change the specifications under SplitIntervals but for whatever reason I keep getting the following error
***********************************************************************

A USER ERROR has occurred: but no positional argument is defined for this tool.

***********************************************************************

Why is this error occurring?
The documentation I am using is here: https://gatk.broadinstitute.org/hc/en-us/articles/360036899592-SplitIntervals

What are all the \% in your output file name? Does it work if you try a simple output file name (e.g. -O foo.list)?

I can't as it is being used on multiple files.

It is meant to be used to create an interval list with the numbers Interval_list_SampleNameHere_.shard_1ofn_.interval_list

Please answer both questions, not just the last one. And even if it is being used on multiple files, try one so you can debug it. That's the way debugging works: you need to try things. Try changing each of the variables you pass to a simple string and see if the error changes. Try and figure out which part of the command is actually causing the error and for that, you will need to run the command by itself.

Thank you for the help I will try with one file - it does seem to be the -O part of the command I will work on it from here and try with perhaps just one file.

This solution where the -O is just set to the output directory seems the most easy way to handle the issue:
Keep the following variables the same:

singularity_image is the path to the GATK sif file
bam_file is a bam file from a defined bam directory
reference is a path to a .fa file
targets_interval_list is a interval list created by Picard (ending with .interval_list)
NUM_BINS is a manually preset value from 100 to 200.

The code with the modified -O:
  # run SplitIntervals with shard name
  singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \
  SplitIntervals \
  -R ${reference} \
  -L ${targets_interval_list} \
  --scatter-count ${NUM_BINS} \
  -O ${output_dir}

While it doesn't name the shards in a very elegant way it seems to work. The shards are made from the interval list and therefore, should not be named after the .bams that will be processed in other steps that are used (for in this case a GATK g CNV pipeline).


text:
What is the limit of partition incompleteness in Bayesian MCMC?
I am seeking advice on whether or not to incorporate incomplete data into my BEAST analysis.
I have HIV data consisting of two partitions of the pol gene: one that is 1.2 kb in length and present in all isolates, & another that is 2 kb in length but only present in one-third of the isolates. I am wondering if I should include the second partition in my analysis, despite its incomplete nature.
I have found some publications discussing the use of incomplete data in phylogenetic analyses, such as Yuchi Zheng & John J. Wiens (2015), Rasmus Hovmöller et al. (2013), John J. Wiens & Jonathan Tiu (2012). However, I remain uncertain about the best course of action.
If you have experience working with incomplete data in BEAST analyses or any insights on this matter, I would greatly appreciate your advice.
Thank you in advance for your help.
Best regards,

Update
As seen from the answer & comments, there has been a misinterpretation of my partitions' structure & their incompleteness. So hereinbelow is a short description of my data, the subject topic, & a diagram.
There are 1009 isolates of the pol gene of HIV from the Republic of Armenia. They were sequenced from the start of the gene by two means: Sanger & an NGS; which caused two groups of length: 1.2 kb & > 3 kb, respectively.
This way the situation may be demonstrated via the next artificial miniature diagram:

Based on the picture like above, I have thought it rational to divide the data into two staple partitions: the first third of pol & its last part. The first partition has all sequences completed, and the second one has mostly full-length complete sequences for a third of isolates and just NNN…NNNs for the other two-thirds.
Ns are automatically added at the alignment step. Sanger sequences are just short and end from 1200 to 1500. I currently do not have access to the source data. They were provided to me already built and aligned at the reference.
In summary, the question is, "Should only the complete first partition be used in subsequent BEAST analysis or should the incomplete second one be included as well?"

@gringer I have updated the question.

Thanks, that's a great improvement; much clearer now.

The original answer was "criticised" for being obtuse.
The basic very simple answer is:

If the data you present in the diagram (alignment) is left EXACTLY in the format used to make the tree - it honestly does not matter whether you have 1 partition or 2 partitions. It makes no practical difference. The posterior likelihood and stuff will change a bit - basically the result is the same. The reason is given in the very long explanation below. I admit 2-partitions looks cooler and there's reason not to do it, computationally it will be a bit slower.
Partitions I am happy to explain why a partition is used  - but thats different answer or question.

What must not be done with HIV is changing every '-' to '?' or changing every '-' to 'N'. Whether this is 1-partition or 2-partitions. With higher eukaryotes that is fine and that will work very nicely, but HIV no. So Weins et al (2012) can do this (almost certainly have done), but you can't.

I know it seems a very small change replacing every '-' with a '?' or 'N', but what it represents seismic. I have attempted to explain why it is so important below the line, but it's about indels (insertions/deletions). I recognise my explanation may not be perfect.
Confusion The confusing bit is you mention 'NNNNNs' and thats the worry. If NNNNs are changed to '-' thats fine, if they are left as 'NNNNs' that raises the complicated spectre of 'indels' and will get very complicated. It probably can be done but its really complicated.
To make the point clear if I perform a trimal analysis on my alignment
trimal -in infile -out outfile -nogap -fasta

It will remove all sites with '-' in them, i.e. it will remove the entire site, not just the '-' character. I think it's really cool. Anyway, the tree from a trimal alignment will be basically the same as '1-partition' or '2-partitions' with all the '-' retained, even though the alignments look very different. That sounds weird but it is a consequence of how Bayes (i.e. Beast) or ML calculate phylogenetic information, which I tried to explain below. HOWEVER, if there is an alignment with 'NNNNs' instead of '-' that will produce a very different tree.
So for simplicity the answer is that 2-partitions is fine for '-', but not if those '-' are represented as 'NNNNN'. Thus for '-' alignments it looks cool, but makes little difference. 2-partitions can make a HUGE difference with different data particularly with loads and loads of NNNNNs instead of '-'.

Summary HIV is a special case because env in particular - and possibly all the HIV genes - are known to contain extensive indels. Indels in env are so prevalent in HIV they can be approximated via distribution (did Ziheng Yang do that?). So the short answer, Weins et al (2012) can do it because they're dealing with higher eukaryotes and genes with negligible indels (in addition they could remove known indel sites), in HIV you would need independent evidence the missing data site was not going to be an indel (i.e. that particular site has never had a known indel in any HIV data set).
The crux is there is a difference between missing data and absent data. Absent data can be a no-go (below; context depending - some times it makes no difference [immediately below]), it can be detrimental to the phylogenetic signal. Missing data is very interesting.
Basically, if you are going to leave the missing 2kb gene as '-' ... thats perfectly fine. It will not add any phylogenetic information - but it will not hurt the calculation. If 'N' or '?' are used instead to represent a missing 2Kb locus for HIV instead  of '-' thats different.

Weins et al (2012) is a very interesting study.
There are four issues

Sample skews
Absent data which causes blockwise deletions in Bayes and ML
Missing data which might recover phylogenetic information at a given site if correctly annotated.
4.Quality of data

The traditional view is that it is better to remove the site/ or sample, if the quality of data is in doubt or for example an entire gene is missing. The traditional view is it is better to be certain of the data regardless of possible risk of a sample skew. However, sample skews are the Achilles heal of population genetics and a lesser extent to phylogenetics. Thats what Weins et al (2012) is concluding, it's really cool. I like it.
We all know that an analyst can't deliberately cause a sample skew (i.e. you can't omit samples without good reason) - that goes with out saying. Right?
The context of the question is "if a sample has pol and its other gene wasn't sequenced - do you fill in the absent gene with "Ns" or "?" to trap the phylogenetic information in pol and the samples which have data in the gene/locus/partition with limited coverage. It's a really good question. My suspicion for HIV is I'm afraid you probably can't, but that might just be env.
Absent data
Bayesian and ML have no '-' character in their matrix. If an indel is present it will trigger that site to be excluded from carrying any phylogenetic information specifically for Bayes and ML. Thus
                Site1
Sample 1         A
Sample 2         A
Sample 3         A
....
Sample 9999      -
Sample 10000     A

Site 1 will not carry any phylogenetic information for Bayes or ML because the single '-' is not represented in the mutation matrix. The entire site for 9999 nucleotides is excluded. A pairwise analysis would be needed, e.g. distance methods.
Missing data However "Sample 9999, Site 1" could be represented as 'N' instead of a '-' because that means any nucleotide in the matrix. That works! Site 1 is now phylogenetically informative under Bayes and ML. Hooray! However, are we certain Sample 9999, Site 1 is a nucleotide rather than an indel? If indels are rewritten as nucleotides, i.e. an absent nucleotide is converted to a nucleotide of unknown type ('AGCT'), i.e. missing data - non-nucleotide to nucleotide - thats cheating.
This creates issues particularly for HIV env where there are loads of indels. If an indel was replaced with an 'N' - the analysis might work nicely, but there is real risk it is not being true to the data.

Possibly explanation Sample skews in phylogeny are really interesting. What can happen is long internal branches can cause changes to the branching order. However, if those internal branches are shorten by more rigorous taxon sampling then these shortened internal branches can recover the true true. I think it was Jo Felsenstein that discovered this (the guy who invented maximum likelihood as a method for phylogenetics in place of parsimony).
One explanation of Weins et al (2012) is that the 'missing data' samples are shortening the distance between internal branches - by providing a more rigorous sample - and that is resulting in recovery of the correct tree. So it is probably a result inline with Jo Felsensteins original observation.
So if pol never has indels - ever - you can do it. This would also help population dynamics and all shorts of calculations. If there's a risk that you are replacing indels with nucleotides it is very much, much better to avoid it.
Conclusion Everyone in phylogentics knows the issues of sample skews - a possible sample skew is permitted if the data is genuinely never been produced - it's not in existence (obviously). However, it is not permitted if the data exists but the analyst hasn't included it - for whatever reason. Its also preferred if the quality is in question - its better to err on the side of caution.
If you fill in an absent gene with "?" or "N" to obtain phylogenetic information from genes that are present for that sample ... that is permitted providing it doesn't change the possible nature of the underlying data. The rationale for doing this is the risk of sample skews.

So would I adopt a "non-traditional" approach? As a general rule for viruses no and this is a particular special case regarding indels.
If this was a higher eukaryote and like - that critical sample and  that gene - didn't give sequence. This has happened. So would I dump the sample (or in this case the locus) - no I'd caveat the tree and explain and empirically check the assumptions.
If I was in a situation where 2/3 of the data were missing, for a higher eukaryote (DNA with great proof-reading), personally I would do a full scale investigation, i.e. phylogenies with and without the gene/locus/partition in question. That is pretty much was Weins et al (2012) have done.

Thank you very much for your reasoning! My question was mostly dedicated to the whole partitions, not single incomplete sites or indels. My sequences jointly cover only pol, but it also has a bunch of indels. I am intent on not including insertions and exchanging deletions with N.

@M__ I am really thankful for your attention! But please, could you shed light on your final point about indels? Have I understood you correctly that we should exclude all insertion columns and exchange deletions with N (I was intent on doing this initially), shouldn't we? Or we should exclude deletion sites as well, or we cannot analyse HIV fragments susceptible to indels with Bayesian MCMC at all?

@Vovin I've summarised the basic answer above the question. Please check it and if it still not clear I will attempt another explanation.


text:
Is there available REGA or COMET versions for downloading & using on a local machine?
I am currently working on verifying the subgenotyping of our sequences and a subsample from the Los Alamos NL database. My supervisor has recommended using REGA & COMET for this purpose. While our dataset, which consists of slightly > 1,000 sequences, can be easily subgenotyped using the online versions of these tools, the LANL dataset is significantly larger, prompting me to consider local installations. Unfortunately, however, I have been unable to find the source repositories for REGA & COMET.
I did come across a link to the entire REGA DB but I am unsure if it contains the subgenotyping tool or if the last can be installed separately. I am also interested in exploring other locally installable tools of comparable quality.
In my search, I found jpHMM, which is available for download and free use. However, it seems to be less frequently used in publications compared to the aforementioned ones.
I would appreciate any insights on the quality of jpHMM or any other recommended tools for local installation (even as packages, e.g. for Python).
If you have experience with these subgenotyping tools or any relevant advice, please share your thoughts! Your help is greatly appreciated!
Best regards,


text:
How can I make this Biopython program (to correct erroneous barcodes) run faster, and is there any alternative method?
This question has also been asked on Biostars
I am looking forward to getting a valuable suggestion for a bioinformatic problem.
Background:
Currently, I am performing a de novo whole genome assembly. At the stage of barcode correction, I lost nearly half of all the reads due to erroneous barcodes.
Barcodes, in this context, are 18-nucleotide-long DNA sequences that were added to the DNA insert during NGS library preparation. Now, I have a file that contains all the 18-base barcodes used in my projects. These barcodes are important because they are unique for reads that came from a DNA molecule and therefore helpful in generating long-range information.

While library construction, 18-base molecular barcodes were used for the linked-read generation, and the barcodes have a high incidence of "N" in positions between 12 and 18, as you may see below.

There is a very large occurrence of Ns in positions 12 to 18 of the barcodes and for this reason, a large pool of erroneous barcodes were filtered.
Now, I have been asked by the company to do one of two things, to recover the lost reads:

Either replace the barcodes(barcode fastq file from sequencing core) with the correct ones from a pool of barcodes (This is a text file I received from the company that contains the pool of all 18 base barcode sequences)

Since most "N"s are seen in positions 12 to 18 of the 18-base barcode sequences (barcode fastq file from the sequencing core), replace the 11-mers in the barcode file, and replace them with the 11-mers from the correct barcode file (text file).

I have decided to use Python code to execute the first method.
#Define the path where fastq file is stored.
import os
from Bio.SeqIO.QualityIO import FastqGeneralIterator
path=os.path.join(os.path.expanduser("~"), "Desktop", "SVA1_S1_R2_001.fastq")

#create a new file handle to write the file
path_new=os.path.join(os.path.expanduser("~"), "Desktop", "trimmedfastq.fastq")
path_ust_barcode=os.path.join(os.path.expanduser("~"), "Desktop", "UST40MBarcodes.txt")
write_file=open(path_new, "w")

#open the text file containg the correct pool of barcodes

with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

handle=open(path, "r")
output = []

#===================================================
# Return the Hamming distance between string1 and string2.
def hamming_distance(string1, string2):
    # Start with a distance of zero, and count up
    distance = 0
    # Loop over the indices of the string
    L = len(string1)
    for i in range(L):
        # Add 1 to the distance if these two characters are not equal
        if string1[i] != string2[i]:
            distance += 1
    # Return the final count of differences
    return distance
#===================================================

#set correct barcode to "None"
corr_barcode=None
#set the minimum hamming distance to "None"
min_distance=19
#open the fastq file

for title, sequence, qual in FastqGeneralIterator(handle):

    # if a perfect match exists, just output as is
    if sequence in barcodes:
        output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))
        continue

    for read in barcodes:

        h_distance=hamming_distance(sequence, read)

        # since if there was a perfect match we won't be here,
        # if distance is 1, we'll never do better
        if h_distance == 1:
            corr_barcode = read
            break

        if  h_distance <= min_distance:
            min_distance=h_distance
            corr_barcode=read

    output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))

with open(path_new, "a") as wf:
    for line in output:
        wf.write(line)

handle.close()
print("Done")

Problem
I would like to know if this is the only Pythonic way to tackle this problem. This method is dead slow. My erroneous barcode fastq file is 30+ Gb, and the UST10x barcode file (the text file that contains all the correct barcode pool) is of size 750Mb. I have split this text file into 1000 smaller files, each about 750kb and having ~40,000 lines per file. It's been the third day and the program is still running without outputting anything. I am looking forward to getting valuable suggestions here, with respect to the code, especially how I can make the process faster.
NB:  I've tried replacing
with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

with
barcodes = list()
with open(path_ust_barcode, "r") as fh:
   for b in fh:
      barcodes.append(b.rstrip())

Still, the code runs extremely slowly.
I am using an IBM server with a 197GM RAM, and an octa-core processor. In my initial code, it was iterating through all the correct barcodes of the text file (UST10x file), for each erroneous barcode in the fastq file, at the end of every iteration, it opened a text file in append mode and wrote the correct version of fastq reads. Since my server has high memory capacity, I rewrote the code in such a way that it stored all the correct barcodes in the memory, and finally wrote it in a text file, by bringing the write function outside of the loop.
Python is not absolutely necessary. I'm not familiar with C or C++, but that being said, I would really love to get this process done. I'm using an IBM server with Ubuntu (Jammy), and 197GB Ram, and the current local storage has only 150GM free space. if you can provide me suggestions as to how to execute it in C or C++, I can do it.

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.

Try profiling (benchmarking) your code to work out where it is slow.
My other suggestions are:

replace your simple Pure Python hamming_distance function with something optimised like RapidFuzz - see https://github.com/maxbachmann/RapidFuzz

cache the corrections you are making - you have tonnes of RAM so a building a dictionary mapping the fragments with N to the corrected fragment ought to be possible. And then for each new sequence, check if you already solved it and if so reuse the same choice.

Thank you for adding a response. I have posted the same question on Biostar, where you can find the code aligned with your logic. A modified version of the code was suggested by @IanSudbery, and you may take a look at it. https://www.biostars.org/p/9561109/#9561603


text:
Why is there antisense sequence in RNAseq data
I'm looking at RNAseq data from CCLE. The data is paired-end.
Take the cell line Hs578T and the gene HRAS as an example.
The cell line carries a G12D mutation (c.35G>A), so the change in cds is:
ggc ggtgtgggca agagtgcgct g - Wildtype CDS
gAc ggtgtgggca agagtgcgct g - Mutant CDS
 ^

My question is, when I grep the mutant CDS gAcggtgtgggcaagagtgcgctg, I do not get a match in my .bam file. But when I grep the reverse complement (i.e., anti-sense sequence), I get matches coming from both mates.
My question is, why are there reads in the .bam file that correspond to the anti-sense strand of the DNA? Shouldn't all the reads of mRNA be the CDS?

I got the answer -
The punch line is that the sequence in .bam files is ALWAYS '+' strand of the reference, no matter what mate the read came from, or which strand the gene sits on.

Yes, the mRNA results should include CDSs (not only, don't forget you also have UTR sequences which are non-coding so are not part of the CDS), but there is no reason why they should all be from genes on the + strand of the genome. Remember that the distinction between "sense/+" and "antisense/-" is essentially a human invention; it doesn't imply that the cell considers one strand to be in any way superior than the other, and genes are pretty randomly distributed across both strands.
This means that when you align to a human reference genome, you will always have some sequences aligning to the + strand, if the RNAs were transcribed from genes that happen to be on the sense strand, and others will align to the - strand, if they were transcribed from genes found in the antisense strand.
This is completely normal and expected, there is nothing wrong with your data. For example, human tP53 is on the antisense strand (note how the arrows are all pointing to the left):

Can you help me understand this: does the bam file sequence always correspond to the +strand of the reference genome, regardless of which mate it came from?

@Nomad420 yes, exactly. Basically, we arbitrarily decided that "this is the 'sense' strand" and then all coordinates are given with respect to that strand.

Thank you. This clears up a lot of my confusion!

Sequences in a bam ar always in the forward direction, with respect to the reference, even if the original read was in the opposite orientation.  The flag indicates if the read as you see it in the bam is just like the original read, or if the bam is displaying the rev-comp.


text:
How to analyze hashtagged data with scanpy?
I used hashsolo to demultiplex my data and I am a bit confused as to how to continue the analysis.
Following hashsolo my samples look like below:
A1
AnnData object with n_obs × n_vars = 8684 × 18197
    obs: 'most_likely_hypothesis', 'Classification', 'n_genes'
    var: 'gene_ids', 'feature_types', 'genome', 'n_counts', 'n_cells'
    layers: 'counts'

A2
AnnData object with n_obs × n_vars = 9901 × 20862
    obs: 'most_likely_hypothesis', 'Classification'
    var: 'gene_ids', 'feature_types', 'genome', 'n_counts'
    layers: 'counts'

I have a few questions regarding the analysis. My initial plan was to merge the 2 samples. However, when I do that and I print the Classification the merged object seems to add the classification counts because the same hashtag names were used.
For example hashtag 1 has 949 in sample A1 and 1286 in sample A2.
When I do print(adata.obs['Classification'].value_counts()) hashtag1 has 2235 counts. Will that be a problem with the analysis?
I merged the object like this:
# add some metadata
A1.obs['type']="C"
A2.obs['type']="T"

# merge into one object.
adata = A1.concatenate(A2)

adata
AnnData object with n_obs × n_vars = 18585 × 18064
    obs: 'most_likely_hypothesis', 'Classification', 'n_genes', 'type', 'batch'
    var: 'gene_ids', 'feature_types', 'genome', 'n_counts-0', 'n_cells-0', 'n_counts-1'
    layers: 'counts'

print(adata.obs['Classification'].value_counts())
Hashtag_2_Antibody    5134
Hashtag_3_Antibody    3150
Hashtag_4_Antibody    2000
Hashtag_1_Antibody    2235
Doublet                                         1598
Negative                                           9
Name: Classification, dtype: int64

My hashtags also show up wen I do sc.pl.highest_expr_genes(adata, n_top=20, ). Is that a problem?

I would also like to plot a UMAP following leiden for each hashtagged sample but I am a bit confused on how to do that. Normally, I would add information such as sample or type on metadata and color by sample for example. But if I try to color by classification it will merge the hashtags from both samples.
Thank you

You could add something like A1.obs["Classification_obj"] = [c+"_A1" for c in adata.obs["Classification"]] but the classification should also be clear in the merged adata object from the "Classification" and "batch" labels together
More importantly, after running hashsolo, you must delete your hashtag feature_types from adata.var, something like A1 = A1[:, [g for g in A1.var.index if g not in list_of_Ab_hashtags]] . Once you have the cell type classifications and are happy with them (can be worth doing further QC first), you can throw away that Ab count data for the downstream standard GEX analysis

Thank you so much! That really helped. Should I also delete the cells that were classified as doublets and negative?