text:
How can I give different names to files in a directory with a for loop in a bash script?
I'm expecting to get 17 different paired-end fastq files (34 in total), so I want to make a bash script to just run my code through all the fastq files in a directory at once. How can I change the name of the input and output files each time the script runs through each file? So when it moves to the file_002, all names have file_002 at the beginning instead of file_001, and so on. And also, when merging the R1 and R2 reads how can I make that it only merges the correspondant files with a loop? for examples merging only file_001_R1 with file_001_R2, file_002_R1 with file_002_R2, file_003_R1 with file_003_R2, and so on.
This is my code:
for file in directory_name
do
pear -f file_001_R1.fastq.gz -r file_001_R2.fastq.gz -o file_001.fastq
cutadapt -g TGATAACAATTGGAGCAGCCTC...GGATCGACCAAGAACCAGCA -o file_001_barcode.fastq file_001.fastq
cutadapt -g GTGTACAAATAATTGTCAAC...CTGTCTCTTATACACATCTC -o file_001_UMI.fastq file_001.fastq
seqkit concat file_001_barcode.fastq file_001_UMI.fastq > file_001_concatenation.fastq
seqkit rmdup -s file_001_concatenation.fastq -o file_001_unique_pairs.fastq
seqkit subseq -r file_001_unique_pairs.fastq > file_001_unique_barcodes.fasta
bowtie -q --suppress 1,2,4,6,7,8 -x ref_index file_001_unique_barcodes.fasta > file_001_barcodes_allignment.bowtie
sort file_001_barcodes_allignment.bowtie | uniq -c > file_001_barcode_counts.txt
awk 'BEGIN{print "Barcode,TF_variant,Code"}{print $3","$2","$1}' file_001_barcode_counts.txt > file_001_barcode_counts.csv
done

Thank you.

I assume the code is directory_name/*, $file is now by each file within directory_name. At the moment $file is simply directory_name as its written. There'll be bash experts here that will describe the loop for the paired end files.

Hi, please consider upvoting and accepting one of the answers. These are good answers from two different perspectives: one pulling the directory file contents; the second using a numerical loop. Both would work.

Here's what I use - along with an example you can modify as needed:
FASTQ_DIR="/directory/containing/FASTQs/"

for R1_FASTQ in "$FASTQ_DIR"*R1*.fastq.gz
do
                         BASE=$(basename "$R1_FASTQ" _R1.fastq.gz)
                         R2_FASTQ="${R1_FASTQ/R1/R2}"
                         R3_FASTQ="${R1_FASTQ/R1/R3}"
                         echo "BASE: $BASE"
                         echo "R1: $R1_FASTQ"
                         echo "R2: $R2_FASTQ"
                         echo "R3: $R3_FASTQ"
                         # process R1/R2 here using "$R1_FASTQ" "$R2_FASTQ" for input, and "${BASE}_XXXXXXX.fastq" for output
done

Basically, the strategy is to iterate over all the R1 FASTQs in a directory, calculate a BASE name that can be used to rename output files with arbitrary extensions (like .sam, .bam, etc...), and use BASH shell parameter expansion to calculate and set variables to the R2 (and R3 if necessary) reads present in an Illumina run that uses standard Illumina naming conventions.
This may need modification depending on the naming conventions you are using.

Something like this? This uses printf to create a zero-padded 3-digit number corresponding to 17 numbers from 1 to 17:
for fileBase in $(printf 'file_%03d\n' $(seq 1 17));
   do echo ${fileBase}_R1.fq.gz ${fileBase}_R2.fq.gz;
   pear -f ${fileBase}_R1.fastq.gz -r ${fileBase}_R2.fastq.gz -o ${fileBase}.fastq
   ...
done


text:
Gene Expression Database for Psychological Disorders
Is there a comprehensive database available online of gene expression across the brain taken from individuals with psychological conditions? Resources like the Allen Institute's Human Brain Map (link here for reference) feature heatmaps of gene expression taken from healthy volunteers, and I'm looking for similar datasets from patients with particular conditions, especially sensory processing disorders such as ASD and ADHD. Any and all data is appreciated!

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

Welcome to the site. To be honest I would suggest you are likely best placed to answer this question, given your expertise demonstrated on other SE sites. NCBI SRA will host the datasets required and their reference number will be given in the associated publications. This is a standard approach to many data sets associated with RNAseq data, i.e. they are exclusively on the established generic repos.

source: Stack Exchange

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

sort: votes


text:
Filtering criteria for non-coding features with very low counts
I am trying to do DE analysis of non-coding features of A. thaliana. I find in the miRNA and lncRNA counts file that they are abundant in zero counts, and most of the non-zero counts are very low. Now, for the DE analysis of coding features, an inbuilt filtering method to remove the extremely low counts is available in edgeR, called filterByExpr(). Is this function suitable for only those features with high counts like coding sequences, or is it suitable for non-coding features with low counts and low library sizes too?
If not, then what kind of filtering criteria should I use for filtering out genes with insignificant counts? Or should I filter out only the ones with row sums 0?

Partially answered before. I am assuming this is a lncDIFF calculation. For the lncDIFF paper it is unlikely to work as the following figure displays,
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-5926-4/figures/4
However, don't take my word for this, test it. You should test the option explicitly via glmQLFTest and glmLRT. If you perform either of these tests with or without the filter the better fit (one not rejected) is the correct answer. I strongly predict the filtered data will be rejected and unfiltered data  (hopefully) will pass.
The question of what is better glmQLFTest versus glmLRT is discussed in the previous answer, here Statistical methods suitable for DE analysis of non coding RNAs
Just keep in mind this is plant data and these tests are developed for humans.

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


text:
DE analysis tool and method for non-coding features
I am currently working with the non-coding features of A. thaliana, and trying to get the DE features. Among the three DE test methods in edgeR such as exactTest(), glmFit() and glmLRT(), glmQLFit() and glmQLFTest(), which is/are the one(s) most suitable for running DE tests on miRNA, snRNA, snoRNA, pseuogenes, transposons and 5' and 3'-UTRs?
Should I be using any other DE analysis tool instead of edgeR for some specific non-coding feature, like there's lncDIFF for lncRNA?

Regarding what are the differences between these tests, this answer by one of the authors might be helpful: https://support.bioconductor.org/p/76790/

The post here is useful ... https://support.bioconductor.org/p/76790/
ExactTest is really cool and completely different from glmQLFit glmFit etc ... you'll need two controls to check it works (it could just reject everything). Its described here https://rdrr.io/bioc/edgeR/man/exactTest.html
There is a concern, some of the lncDIFF data looked more like a lognormal to me in the publication, but I didn't read it in sufficient detail. Anyway your data will probably approximate to a negative binomial distribution and it's statistically tricky to  assess whether two negative binomial distributions are significantly different ... the ExactTest will do this.
Is the total output of the control significantly different to the test sample? ... the ExactTest will tell you. You want p < 0.05, preferably a lot lower than 0.05 to be convinced. However, you need to compare two control samples against each other, where p > 0.05 for the subsequent tests to be valid.

So really nice BUT you need to check whether it is a negative binomial FIRST. Thats what the other tests are doing. Please do check whether lncDIFF identified a negative binomial, because if it ain't ExactTest will not work, they could be using a different distribution.
The other tests you mentioned are different to the ExactTest because they simply assess whether the data is negative binomial. The distribution is the frequency of each RNA transcript BTW. edgeR's lnfDIFF requires that distribution and thats how it can assign probabilities to the difference between the frequency distribution of two identical transcripts under difference experimental conditions (or other conditions). So you need a test to say "yeah I'm a negative binomial" (assuming it is negative binomial), then everything else via lncDIFF is then valid.
First you've got to fit you data to a model ...
y <-glmGLMFit(...)

... is the function and that uses general linear models to fit the negative binomial. The fit will result in something called "a residual" and that is very important, the smaller the better. Testing whether that residual is small enough for a negative binomial to be declared is down to ...
glmLRT(y ...)

This is a likelihood ratio test, to test the observed likelihood against the model (p > 0.05). You can also use
gof(y ...)

this is goodness of fit, as the link says the residual follows a chi-squared distribution, i.e. its tested on the chi2 tables and thats easy to understand. Again you need a probability > 0.05 to declare a fit.
HOWEVER, the lnfDIFF guys developed a quasi-likelihood to assess the distribution
z<-glmQLFit(...)

Fits the data via quasi-likelihood, I honestly have no empirical idea how that works.
glmQLTest(z ...)

tests the fit to the distribution. You need a probability > 0.05. Again I have not empirical idea whats going on here.
ANY non-significant result is good, doesn't matter if its GLM with gof, GLM with LRT, GLM-QL (whatever that is) because all of these are solid tests.  Which would I use? The guys who developed lncDIFF used QL (quasi-likelihood) to assess the fit (might be a good idea). However, statistically I honestly don't know whats going on, so I would also look at the GLM and fit it against the goodness of fit - because thats easy to understand.
What do I do to assess the fit? I use machine learning, but it's way more complicated than anything here.

From the comments.
It's important to note the context. The OP is performing non-coding RNA work, lncDIFF is much better than normal edgeR for this type of data described here Statistical methods suitable for DE analysis of non coding RNAs.
lncDIFF is now available for R here
... so the methods of edgeR should be available to it. glmQLFit() is a general function in R. Providing the lncDIFF object will be accepted by the tests, by performing everything in R this should solve this.
I would advise using lnfDIFF's QL directly, i.e.
y<-ZIQMLfit(...) # fit the distribution
ZIQML.LRT(y ...) # test the fit

The manual is here
https://github.com/qianli10000/lncDIFF/blob/master/lncDIFF-manual.pdf

I hope that makes sense these are pretty advanced stats so it can be tricky to understand them. Whilst I can teach code, I don't teach complex stats outside machine learning.

So... is edgeR a good tool for analysing the DE of noncoding features? I mean since lncRNA has a separate specialised tool for its analysis, which they claim is better than edgeR, is edgeR fine if I want to analyse noncoding features?


text:
is it possible to count cell in Violin plot in seurat?

I tried to quantify the number of cells of Fzd9+ Macrophage.
how can I know the cell count?
thank you

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

Does @haci's response answer your question? If so please mark it as "accepted". The question appears to contain very little detail and is easy for a misunderstanding to occur in these situations.

This can be a good starting point: table(your_seurat_object@meta.data$the_column_that_has_the_labels_of_interest)
Basically you will be counting the number of rows (each single cell is a row) in the meta.data slot of your Seurat object for each level of your column of interest (for example cell type annotations).


text:
BioPython SeqIO fails with Invalid Alphabet found
I am trying to store a GenBank file that could be read by Geneious later.
I found an example how to create the GenBank file here
But now I run into Alphabet issues.
I tried all protein classes from the documentation, but all fail with the same error:
TypeError: Invalid alphabet found, <class 'Bio.Alphabet.ProteinAlphabet'>.
I dont know what else to do here? Using biopython==1.73
    record = SeqRecord.SeqRecord(
        Seq.Seq("MIRQALAVAALLLAGTAQADGLIDN", alphabet=Alphabet.ProteinAlphabet),
        id="FooBar",
    )

    with NamedTemporaryFile(mode="w") as genbank:
        SeqIO.write(record, genbank, "genbank")
        pass

If possible, please show a minimum reproducible chunk of code. The code shown here does not include any library imports, so will not work as-is.

Quickly ...
from Bio.Alphabet import generic_protein
from Bio.Seq import Seq

then ..
myprot = Seq('EXAMPLEPROT', alphabet=generic_protein)

Thanks, I tried without setting the alphabet but Seq needs that parameter. Let me try another way of importing the Alphabet

Using generic_alphabet tells me this: ValueError: Need a Nucleotide or Protein alphabet

I was wrong, generic_protein does the trick, not generic_alphabet

Thanks @ElDude done .. many thanks.


text:
All against all sequence alignment
I have a large set of protein sequences and would like to create a "similarity matrix" out of them (for phylo trees later). So I have to run n(n-1) sequence alignments. Wondering what library I should use for that and also what sophisticated lib would provide this all against all without me writing a for loop and parsing results myself...
Currently using Bio.pairwise2.align.globalxx from BioPython as a development placeholder but I do not think the results look really useful? It returns a score > 1 and no similarity in [0,1[ or evalue.
I have been out of the Bioinfo game for a while, so any suggestions for libraries (python preferred) much appreciated! :)

I wouldn't use Python here. I've looked at using Bio.pairwise2.align. and its has major shortcomings. In this context the problem is establishing homologous alignment positions ... if you want to make a tree.
Whats the problems? This is you start introducing the gap character '-' as an additional state in constructing your tree and the tree will bias towards not the maximum number of substitutions but the maximum gap (indel) difference between the sequences.
The second problem is changes in homologous sites between different pairwise alignments. This gets complicated fast.
Is that bad? The problem is we have little understanding of the mutation model behind generating indels and there isn't a model to account for varying descriptions of homologous alignment positions. However, we have a very clear understanding of amino acid substitutions in a single global alignment, i.e. a singular description of homologous amino acid positions. A lot of work has gone into modelling amino acid mutation, with a tree a practical output of that work.

Thus, I would simply align using muscle5 under the -super5 option. This will produce a global alignment very quickly and for vast amounts of data. This alignment will produce the best model of homologous amino acid sites.
If you just want a distance matrix then import the entire alignment into MEGAXI and request the JTT model of amino acid substitution with the gamma distribution. This will output a distance matrix nicely and will then allow you to produce a tree.
BTW estimating the gamma distribution without maximum likelihood can be difficult in distance methods (MEGA might have resolved this). This tree will be much better than the one you are proposing (you can estimate this parameter via e.g. IQTREE or ModelTest).
Why? Again alot of work has been done on modelling amino acid substations and JTT is one of the better models (assuming this is nuclear sequence not mtDNA). The LG matrix is much better BTW, but I don't think thats in MEGA.


text:
Metrics to use to estimate variant library evenness, uniformity, or bias
I have a 15k library of variants (OTUs) from a genome region. I want to compare the evenness of this library and understand of the library bias.
I tried different methods to calculate evenness or uniformity of this library such as Shannon evenness index, relative abundance rank plot, kurtosis , skewness. I am confused with the results and why they are not aligned. For example my Shannon evenness index is 0.98 which means library is very even, while my kurtosis is very large (120), indication not uniform distribution.
I should note that for calculating kurtosis I use relative abundance which is read count of every variant /total read count of all variants in that library. My question is what metrics I can use to have a better understanding of library evenness or uniformity, other that above metrics, and why my results are so contradictory?

Shannon and Kurtosis are very different statistics. Kurtosis is generally used in conjunction with skewness as a very general approach to normal distribution signal. I would simply say they are measuring different things


text:
GEO Methylation dataset "GSE73303" shows 0 features
I'm working for a project using GEO, in order to analyze some gene expression data. The SuperSeries I am referring to is GSE73304, but I found some problems working on the SubSeries GSE73303.
I started my analysis Microarray expression dataset and everything was fine. Now I have to investigate the other dataset, the one target as Bisulfite-Seq. But when I try to see inside it, I find there are correctly 20 patients, but 0 features, so seems it is empty. I attach R output:
library(GEOquery)
gsemet <- getGEO("GSE73303")

> gsemet
$GSE73303_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 0 features, 20 samples
element names: exprs
protocolData: none
phenoData
sampleNames: GSM1890515 GSM1890516 ... GSM1890534 (20 total)
varLabels: title geo_accession ... Sex:ch1 (44 total)
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
pubMedIds: 27437034
28883895
Annotation: GPL11154

Could someone help me to find where I'm wrong? Or maybe, if the file is empty, pass me the right file?
Thank you very much in advance

It's NGS data, while GEOquery only returns data for arrays. You will either need to use the data they provude in the supplement of the GEO submission (if they do, it's not mandatory), or download and process fastq from scratch.
It's not gene expression by the way, it's methylation data, bisulfite sequencing.


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

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

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

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


text:
Read methylation data from tar files on GEO (GSE73303)
I'm trying to analyise methylation data using GEO, in order to conduct a multi-omics analysis. The SubSeries I am referring to is GSE73303. Unfortunately, the code used to inspect it produces empty file, so I downloaded tar files, as I found on bioconductor forum (ref. this).
I was able to download the tar files and extract thenm but afterwards I got a list of .bigwig files. I tried to read them on R, but actually I feel very lost and I can't understand how to trace back methylation data from these files.
If someone can guide me into this process, it would be really appreciated
Thanks in advance

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

Ok I modify it based on your suggestion, thanks

As I've already told you here GEO Methylation dataset "GSE73303" shows 0 features GEO/GEOquery does not provide ready-to-use processed NGS data. It can be that the authors provide processed data, but if not (apparently not the case in this submission) then you have to start from the fastq files. Please see the many tutorials online on downloading fastq and processing bisulfite data, or try to use existing pipelines (nf-core methylseq might be an option) or seek a collaboration.

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

@Riesz98 it would be etiquette to upvote.


text:
What kind of data do I need to be able to complete Differential Expression using Limma
I am very new to bioinformatics but have some background in coding with R. I have been asked to find a data set and complete differential expression using the Limma package in R.
I have tried a few data sets but haven't been able to get any code to work.
What kind of data do I need to have in order for Differential Expression to work with Limma R? I have tried load and read.table, is there a read function that I am missing?
I was advised to look on Gene Expression Omnibus to find published data sets to work on but again I am struggling to find any that will work.
Thank you

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

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

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

Limma is relatively generic in what it accepts. Typically it is either data on log2 scale that are either already normalized or need some of the inbuilt normalization methods (such as quantile normalization), or in case of RNA-seq you can provide raw counts which then can run through the limma-voom pipeline.
The format is always a matrix of counts/values, with genes/observations in rows and samples/specimen in columns.
For practice, you could follow the RNA-seq/limma-voom section (chapter 15) of the limma manual which has a RNA-seq case study available at https://subread.sourceforge.net/RNAseqCaseStudy.html with example data and code.
GEO has a lot of published data, often microarray gene expression values or RNA-seq, but lets discuss this in a separate question if you need it to keep this here focused on limma.

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


text:
Sql query slow on Chembl postgres database
Doing complex sql query on Chembl database and it takes 40 minutes on the machine from lab.
I was thinking using dask, extract each table and export to parquet, when do the data processing .
Is this a good approach or there is a better alternative?

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

I suspect the issue is you want to parallelise (multi-threading really) the job to reduce run time.
parquet is recommended. It is simply an archive format that allows compression and is often used with a pandas dataframe. Use parquet to store data by all means, but thats easy via pyarrow and the pandas command
    df.to_parquet('/User/OP/location/file,parquet' ,engine='pyarrow', compression='gzip')

You will need to install pyarrow via pip or conda, there is a faster version (forget its name, something like fastparquet or something used by default). Yeah it's a very good archiving format, in fact its industrial standard and I should use it more.
Python has loads of ways to parallelise of which dask is one way. dask is an entire framework. To leverage dask to parallelise postgres I would refer to this post.
I don't use dask, I am aware it is a solution in high performance Python. Personally I would use this approach directly via multiprocessor
from psycopg2 import connect
import multiprocessor

... etc ...
I would suggest its horses for courses, I simply know how to multiprocess in Python so I don't bother with dask.
To answer the question in summary

parquet is the best, if SPEED is not important. If speed is important then feather
dask is mainstream within high performance Python and its doable. There's no reason not to
There are other multithread Postgress approaches and thats what I would do, but if you just getting into parallelisation you may as well use dask. It will take some learning though.

If you're not into parallelisation and you're not writing a data pipeline for publication I would simply consider getting a machine with more RAM. It shouldn't be taking that length of time, parallelisation will solve the run time - no question - but it's the time it will take to get to grips with a parallel solution. What is likely occurring is the job and database size is exceeding the RAM capacity of the machine. At a wild random guess a 16 or 32 Meg RAM chip would collapse the runtime to a few minutes and thats what you are looking for.

From the comments.
parquet is really for data compression and archiving. The issue is I/o is reasonably fast but not lightening fast. It is probably slower than pickle (it will be slower if you're compressing via gzip as well).
For a 'temporary' save (couple months) and really fast i/o then its feather. There is no dispute about that.
If you are keeping your data a long time and need to save space its parquet. Performance-wise feather is way faster than parquet. I've never used feather BTW I use pickle or parquet
To give some idea a 500Meg flat file compressed an archived via parquet resulting in a 50Meg file, and will take several minutes to write/import. Thats 'cause it's doing loads of stuff. My suspicion is that this too slow for you.

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

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


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

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

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


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

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


text:
how to install cibersort in R, to do deconvolution using my RNA-seq Data with public single cell data
I have RNA-seq Data from lung cancer immunotherapy patients but i don't have any single cell data from my sample, so i wanted to do a deconvolution with a public single cell data. I thought of doing the analysis with the cibersort website plataform but i wanted to obtain the results in R, to be easy to work with them.So i followed the steps in this article, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5895181/ :
library(e1071)
library(parallel)
library(preprocessCore)
source("CIBERSORT.R")
results <- CIBERSORT(‘sig_matrix_file.txt’, ‘mixture_file.txt’, perm=100, QN=TRUE)

But in when i did source("CIBERSORT.R"), i obtain this:
"Error in file(filename, "r", encoding = encoding) :
  cannot open the connection
In addition: Warning message:
In file(filename, "r", encoding = encoding) :
  cannot open file 'CIBERSORT.R': No such file or directory"

I don't understand i was supposed to have a file?

Is the R script available on the Cibersort website? I was looking for it some time ago but couldn't find. I don't have my credentials for the site on this computer so can't check now.


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

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

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

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

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

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

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


text:
Enrichment analysis of DMPs | new Ilumina EPIC V2 900K
Im analyzing methylation data from the new illumina's EPIC V2 array (900K).
I would like to perform enrichment analysis on DMPs. In the past, I have always used missMethyl but since there are a lot of differences between EPIC V1 and EPIC V2 arrays I'm unsure if missMethyl can be used (Im afraid results might be biased)
Does anyone know any enrichment tool already available for the EPIC V2 array or alternatives to missMethyl?
Thanks

We have been using SeSAMe, it works well for the most part. They have an enrichment tool called KnowYourCG. You can check it out their vignette: https://www.bioconductor.org/packages/release/bioc/vignettes/sesame/inst/doc/KYCG.html.
Note, SeSAMe and Minfi (and that suite of packages) have different data structures; However, for enrichment analysis, you only need probe IDs. Then you can use either curated databases from SeSAMe, or you can simply extract gene information, and run gprofiler2 as shown at the bottom part of the link I shared. Hope this helps.

Welcome to the site and thanks for your answer


text:
How to import the average expression matrix for all detected genes in each cluster into Seurat？
I download a matrix which has 19 columns and the first column is gene names, other is cell cluster:
In fact, I just want to plot the feature plot which like this:
My idea is to use the function DotPlot in Seurat package to plot it. So how can I import this matrix into Seurat object or is there any other method?


text:
How evaluate biclustering algorithm over real datasets containing labels?
I have a set of real data containing labels. However, in some biclustering algorithms (e.g., Cheng and Church Algorithm), originally authors applied gene expression data like Yeast without having labels. My question is that what are appropriate evaluation measures in this scenario?

Welcome to the site. I understand machine learning (ML) but the question is not clear. The key Issue is the use of links. We will not click two  external links to understand the question, so it needs a summary of the intended algorithm is needed within the question. Biclustering sounds to me to be unsupervised learning and that would explain the difference. More generally, explaining what data you have and what the labels are is usually helpful.

Biclustering or co-clustering is an unsupervised learning task with which both rows and columns in a data matrix can be placed in concurrent clusters whereas non-overlapping clusters are generated in clustering algorithms.
Data set - a resource of regulatory elements (REMs) contains features like Gene ID, REMID, Genomic location, Predicted function (acts as a label having values such as activation, repression), Model score (normalized absolute binary logarithm of the P-value).

Briefly evaluation measures, which I presume means "features", are determined using supervised learning using a technique called "feature selection". Unsupervised learning generally is used to identify categories for a supervised learning algorithm. The orthodoxy has changed.
"Unsupervised learning feature selection"

To identify suitable features using unsupervised learning a priori then genome location which must have a major outcome, thereafter it would be a case of deleting a given feature and re-evaluating the clusters. If the clusters disintegrate then that is a key feature (not very cool but would work for datasets with limited numbers of features).
Weighted unsupervised learning is popular, this is basically "feature selection", i.e. which features contribute to the clustering. This is basically the formal version of point 1. Personally I have had little success with this approach, i.e. the results were not informative, but it does work. The results are not discriminatory between different conditions, which is what I was seeking.

Regarding labels, labels are not used in unsupervised learning approach - except weighted unsupervised learning - hence their presence or absence usually makes little difference. The only "label" used is the sample type, i.e. the "id" and not the "feature". Without intrinsically understanding the experiment it is difficult to forward an accurate a priori assessment.


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

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

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

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


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

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

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

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


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

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

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

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

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


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

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

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

The graph created from this is the following:

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

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


text:
WGCNA soft threshold selection
Here I am conducting WGCNA analysis. When calculating the soft threshold, it is basically below 0.8. I choose Power=6, which only has 0 and 7. Is it feasible?

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

Please copy-paste text as text into question posts. This makes the text easier to read on some devices, and allows it to be more easily modified when required.

I don't know the context (how many samples, what tissue, experimental design etc) but the mean and median connectivities are pretty high at power 6. I would first check the data for strong drivers of correlation and possibly adjust for them if they do not seem biologically relevant to the question you study. Assuming you find no issues, I would go with a higher power here, maybe 12.


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

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

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

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

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

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

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


text:
Is a classification tree appropriate method to use for my project?
I have about ~3000 genomes of organisms and presence (recorded as 1)/absence (recorded as 0) data of ~300 genes in them. Among 3000 genomes, we have experimental data on 20 organisms, which were represented by 2 different phenotype groups (a, b). According to a previous study, the 300 genes I mentioned are believed to be associated with phenotypes 'a' and 'b'.
To predict the phenotypes of 2980 genomes with no phenotype information, I did hierarchical clustering of 3000 genomes based on the presence/absence of the 300 genes. I cut the resulting dendrogram at a certain height that allowed distinction of clusters that included either group 'a' or 'b', but not both. This gave me 120 clusters, including clusters with no group 'a' or 'b' organisms, clusters with at least one group 'a' organisms, and clusters with at least one group 'b' organisms. Among these 120 clusters, there are 15 clusters that had either group 'a' or 'b' organisms.
I am questioning which collections of genes are responsible for causing phenotypes 'a' or 'b'. To address this, I am thinking about assigning different names for the 15 clusters that contain group 'a' or 'b' organisms, and make a classification tree using these 15 clusters. Does this sound like an appropriate application of classification tree approach? If not, could there be a more sound way to address this question?
Thanks!

This isn't a cool strategy because dendrograms are very poor representations of the relationships between genomes, for example they use UPGMA style metrics fail to account for mutation rate variation between genomes and prone to long-branch attraction artefacts.
I would recommend doing whole genome alignments using e.g. muscle5 -super option and then performing a formal phylogenetic tree using e.g. iqtree or even just MEGAX (distance methods) and using bootstrapping to identify the robustness of clades comprising a and b.
If two clades are supported by bootstraps > 80% then your hypothesis looks fine. At present your downstream analysis could be an artefact of the dendrogram phylogeny.


text:
How to classify TF motifs by family
I have a FIMO output, the input sequences were putative promoter sequences of several species. I want to graph the positions of these motifs along a horizontal graph, but I want to filter the data by clustering the motifs by family, e.g. all the zinc finger motifs being represented by one color, all the helix loop helix motifs being represented by one color, etc. How would I go about clustering the raw .tsv file?

Please explain what FIMO is, and show the first few lines from your raw .tsv file. What language / tool would you like to use for clustering? Is this something you'd do directly in the FIMO tool, or externally using another program?

JohnDoe23, are you looking for an output like the one shown here: https://bioinformatics.stackexchange.com/a/20132/3967

How do you know the "family" to which your motifs belong?  If the source for the motifs that FIMO matched is, say, JASPAR Vertebrata CORE, then you might reasonably decide to use JASPAR's "structural class" as provided by their pre-computed "Annotation table" in JASPAR CORE Vertebrates clustering.
Assuming you are using JASPAR and you do make this decision...
Your "raw .tsv file" is in FIMO Results TSV Format.  "Joining" it with a downloaded copy of the "Annotation table" serves to annotate it with additional columns, having the effect you seek of clustering the FIMO results on the "structural class" of each matched motif.  Joining can be implemented following Running queries directly against CSV or JSON using the sqlite-utils command-line utility.
How you would then plot these features and encode the class as a color requires more detailed description what the plot is supposed to represent or look like.  For instance, do you "want to graph the positions of these motifs along a horizontal graph" with coordinates taken with respect to the TSS of the transcript from which the promoter is taken?

Comment: I don't understand what you mean by "I want to graph the positions of these motifs along a horizontal graph"? Are you talking about a Manhatten plot? Or something else? And I'm not sure about "filter the data by clustering the motifs by family". What type of clustering do you want to perform?
Without further information/description it's difficult to understand what you're trying to do. If possible, please edit your question with an example of what your desired output. Also please indicate the language you're using for your visualisation, e.g. R, python, other.

Answer: I may have completely misunderstood, but perhaps this approach could be adapted to suit. Using the example FIMO output file from https://meme-suite.org/meme/doc/examples/fimo_example_output_files/fimo.html?man_type=web:
library(tidyverse)
library(qqman)
#>
#> For example usage please run: vignette('qqman')
#>
#> Citation appreciated but not required:
#> Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731
#>

fimo_df <- read_tsv("~/Desktop/fimo.tsv", comment = "#") %>%
  rename("SNP" = motif_id, "CHR" = sequence_name, "BP" = start,
         "P" = `p-value`, "zscore" = `q-value`) %>%
  mutate(CHR = case_when(str_detect(CHR, "X") ~ "23",
                         str_detect(CHR, "Y") ~ "24",
                         TRUE ~ CHR)) %>%
  mutate(CHR = parse_number(CHR)) %>%
  select(-c(stop, strand, matched_sequence))
#> Rows: 1494 Columns: 10
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr (5): motif_id, motif_alt_id, sequence_name, strand, matched_sequence
#> dbl (5): start, stop, score, p-value, q-value
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

manhattan(fimo_df, p = "zscore", logp = TRUE, ylab = "q-score",
          genomewideline = FALSE, suggestiveline = FALSE,
          main = "Manhattan plot of TF motifs",
          chrlabs = c(1:19, "X", "Y"))

If you want to colour the points by TF motif (or by 'cluster'), perhaps ggplot would be appropriate, e.g.
library(gghighlight)
axis_set <- fimo_df |>
  group_by(CHR) |>
  summarize(center = mean(BP, na.rm = TRUE))

ylim <- fimo_df |>
  filter(P == min(P)) |>
  mutate(ylim = abs(floor(log10(P))) + 2) |>
  pull(ylim)

ggplot(fimo_df, aes(x = BP, y = -log10(P),
                    color = as.factor(CHR),
                    size = -log10(P))) +
  geom_point(alpha = 0.75) +
  # scale_x_continuous(label = sort(axis_set$CHR), breaks = sort(axis_set$center)) +
  scale_y_continuous(expand = c(0, 0), limits = c(3.95, ylim)) +
  scale_color_manual(values = rep(
    c("#276FBF", "#183059"),
    unique(length(axis_set$CHR))
  )) +
  scale_size_continuous(range = c(0.5, 3)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
  ) +
  gghighlight(motif_alt_id == "TBP")
#> Warning: Tried to calculate with group_by(), but the calculation failed.
#> Falling back to ungrouped filter operation...

Created on 2024-02-05 with reprex v2.1.0
Does that help at all?


text:
How to visualize FIMO outputs?
I used FIMO to scan several of my promoter sequences for TF binding motifs. The output .tsv file gives me the motif name/sequence plus the start and stop location, how can I visualize these outputs to get a graph that shows each motif as differing length rectangles along each promoter sequence, along with a gradient that shows the p-value?

If the sequences you are scanning with FIMO are supported by the UCSC  or Ensembl genome browsers you can upload the GFF file generated by FIMO as a custom track: UCSC Custom Track, ENSEMBL Custom Track. Note that for the Ensembl website the FIMO GFF file is in GFF3 format.

Potential workaround: use MAST.
As far as I understand, 'under the hood', so to speak, FIMO and MAST work pretty much the same. The difference might be just that FIMO focuses on scoring motifs, while MAST is scoring sequences (https://groups.google.com/g/meme-suite/c/qMSj6AdFxoM)
Example of FIMO output: https://web.mit.edu/meme_v4.11.4/share/doc/examples/fimo_example_output_files/fimo.html
Example of MAST output: https://web.mit.edu/meme_v4.11.4/share/doc/examples/mast_example_output_files/mast.html
There might be small differences in the output of the tools, but then it's up to you to see whether this is relevant for your data.

FIMO and MAST have entirely different purposes. FIMO identifies statistically significant matches to motifs in a collection of sequences. MAST identifies the sequences that have the best overall match to a collection of motifs. The output of MAST does include a cartoon showing the motif matches for a particular sequence, it only does this for the sequences that have the best overall match to ALL the motifs in the input. It won't necessarily show the  best matches to individual motifs. The typical use case for MAST is identifying promoter sequences regulated by multiple transcription factors.

(cont) While internally MAST does perform a step equivalent to FIMO the final MAST output is not equivalent to the FIMO output.


text:
Identifying the promoter on the reverse strand using visualization software (seqmonk)
I feel like this might be a really silly question and I've been going back and forth for awhile and I think confusing myself- so sorry if this is dumb.
I am using Seqmonk to visualize genes in an assembly (but I suspect the issue would be the same in any visualization software). I am trying to identify promoters of specific genes to determine if they are differentially methylated. Genes on the forward strand (highlighted in red), the promoter would inherently be upstream- in this case to the left (i.e. if my start site was 5000 and the end site 6000, I would assume I could find my promoter to the left of the start site between 3000 and 5000 bases.
My confusion is what happens with genes that are on the reverse strand (in blue here). Start site and end site could be the same (5000 and 6000- as end sites are always bigger than the start site). Does this mean that in the software, the promoter is expected to be on the right hand side (i.e. bases 6000-8000)? This has been my going thought but now I wonder if having bases that seem downstream (because they're bigger) makes sense, even though it really does because downstream on the forward strand would be upstream on the reverse strand. Or if I need to always do 3000-5000 bases as the promoter, regardless of the the strand of origin.

If I'm looking at a gene on the reverse strand, where would I find the promoter.

If you intended to post an image, we are not able to see it. Not clear what red and blue are.

Have a look at this UCSC browser shot.

Looking at the WDR92 transcript, note that it is transcribed right to left. This would indicate that the promoter is on the right. Indeed, we see a variety of promoter histone marks on the right of the gene, but not on the left.
In other words, yes, for genes on the reverse strand, the larger coordinate will be at the beginning of the gene, and the promoter will include a larger coordinate than the rightmost transcribed base.

THANK YOU. This is what I thought but I suddenly got wrapped up in my head. Super helpful.

@user18170 No problem. consider upvoting/accepting the answer if you found it helpful, in case others find this in the future.


text:
I am trying to create a subset of 10k variants from 25-30 unmapped contigs of a g.vcf file including the header
My objective is to take a g.vcf.gz file and from 25-30 unmapped contigs with titles like "NW_020192317.1", I want to make a subset of ~10k variants from each of the unmapped contigs and make one final g.vcf file that includes the header from the original g.vcf.gz file.
From the post here, I have tried running:
samtools view -bo output.g.vcf -s 123.4 inputfile.g.vcf.gz NW_020192317.1

and I keep receiving the output:
Aborted (core dumped)

I'm not too picky about the number of variants per contig, as long as it is in the thousands, that's why I tried using ".4"
Alternatively, I tried using the code:
samtools view input.g.vcf.gz "NW_020192317.1:1-2" > output.g.vcf

and I recieved the same response of "Aborted (core dumped)" with nothing else. I am confused what is wrong here and why it is throwing an error.

You at least need to be using bcftools, not samtools, if you are working with vcf files. samtools is designed for bam/sam, not vcf, and that’s why it’s throwing that error

Do you have a list of contigs of interest? Are those contig names the first field of the VCF? Does this need to be a random selection or can you keep the first 10k variants from each contig?

Also, gvcf files contain both variants non-variant lines. Do you only want variants?

After some guidance from my advisor, this command allowed me to get the output file I needed.
bcftools view -o newvcfFileName.g.vcf -s 123.2 locationOfTheFile chrom1 chrom2 ... chrom23(mito DNA) --force-samples

I listed every chromosome's title until I typed --force-samples. I found it simpler to copy and paste the names in.


text:
How to identify TF binding motifs?
I have DNA sequences from a single gecko species for a set of genes (SULF1, SOX9, SAL1), I wanted to compare the putative promoter sequences for each gene, and compare them to the promoter regions of other related species. Would something like HOMER work for say finding TF binding motifs in around 20 different sequences for 20 different species, or would I need a larger dataset?


text:
Binning then assembling sequencing reads
I am attempting to do something similar to the "Phylomapping"  pipeline described in this paper: https://www.nature.com/articles/s41559-019-0914-2. I have sets of Illumina PE reads, that I want to compare against a reference genome that is annotated for CDS and CNE. So far i've used makeblastdb to create a custom db for the reference genome, and i've blasted the forward and reverse reads that I have against the reference genome with blastn and dc-megablast. I'm not sure the approach I should take to actually sort or bin the reads, and then subsequently assemble with CAP3.

Please explain what you've already tried. The figure you've displayed mentions contig assembly with CAP3 and UCLUST, and binning with NextGenMap. Have you tried using these programs? If so, what specific issues have you encountered?


text:
how is log2foldchange calculated in DESeq2
I am trying to calculate log2FoldChange in DESeq2 by hand because I am trying to reproduce my results and understand how DESeq2 works.
What I want to know is:

Are the counts first fitted to a negative binomial? If so, how can I get access to those values?
Regarding log2FoldChanges, I read that the coefficents are the estimates for the log2 foldchanges for each sample group. However, these estimates do not account for the large dispersion we observe with low read counts. To avoid this, the log2 fold changes calculated by the model need to be adjusted.
Why? Didn't we just fit the counts to a negative binomial, which should take into account the dispersion.
Finally, how are the log2FoldChanges calculated? It's not possible to figure this out using the raw code because most of the real calculations call C scripts.

Have you checked if the DESeq2 vignette answers your questions? https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

were you able to figure this out? I am trying to do something similar where I have plot showing the differential expression results for a single gene from the DESeq2 output. I plotted the marplot of the log2FoldChange and errors bars with lfcSE. In the barplot I want to include the gene expression ratios that were used in estimating the statistics. Was wondering if it is possible? Somewhat related to your question, but the difference is I am interested in the ratios of the log-counts
Any thoughts or hints will be appreciated.


text:
tensorQTL interaction issue
I hope that the question has its place on this board.
I am trying to run the mapping of cis eQTL using an interaction term with tensorQTL. I am running into the following index error when running the software:
>>> cis.map_nominal(genotype_df, variant_df,
...                 phenotype_df,
...                 phenotype_pos_df,
...                 prefix = prefix,
...                 covariates_df=covariates_df,
...                 interaction_df=interaction_s,
...                 maf_threshold_interaction=0.05,
...                 run_eigenmt=True,
...                 output_dir=output,
...                 write_top=True,
...                 write_stats=True)
cis-QTL mapping: nominal associations for all variant-phenotype pairs
  * 569 samples
  * 24665 phenotypes
  * 68 covariates
  * 69763935 variants
  * including 1 interaction term(s)
    * using 0.05 MAF threshold
  * cis-window: ±1,000,000
  * checking phenotypes: 24665/24665
  * Computing associations
    Mapping chromosome 1
    processing phenotype 137/24665Traceback (most recent call last):
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/core.py", line 197, in calculate_interaction_nominal
    Xinv = torch.matmul(torch.transpose(X_t, 1, 2), X_t).inverse() # ng x (1+2*ni) x (1+2*ni)
torch._C._LinAlgError: linalg.inv: (Batch element 2704): The diagonal element 3 is zero, the inversion could not be completed because the input matrix is singular.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/cis.py", line 294, in map_nominal
    res = calculate_interaction_nominal(genotypes_t, phenotype_t.unsqueeze(0), interaction_t,
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/core.py", line 200, in calculate_interaction_nominal
    i = int(re.findall('For batch (\d+)', str(e))[0])
IndexError: list index out of range

However, when I run the same command but only on 50 genes, the code runs properly.
>>> cis.map_nominal(genotype_df, variant_df,
...                 phenotype_df.iloc[1:50],
...                 phenotype_pos_df.iloc[1:50],
...                 prefix = prefix,
...                 covariates_df=covariates_df,
...                 interaction_df=interaction_s,
...                 maf_threshold_interaction=0.05,
...                 run_eigenmt=True,
...                 output_dir=output,
...                 write_top=True,
...                 write_stats=True)
cis-QTL mapping: nominal associations for all variant-phenotype pairs
  * 569 samples
  * 49 phenotypes
  * 68 covariates
  * 69763935 variants
  * including 1 interaction term(s)
    * using 0.05 MAF threshold
  * cis-window: ±1,000,000
  * checking phenotypes: 49/49
  * Computing associations
    Mapping chromosome 1
    processing phenotype 49/49
    time elapsed: 0.44 min
    * writing output
done.

The interaction factor is a categorical variable (0/1), I am not sure why the issue is being raised and how to fix
I would appreciate any help/direction to go forward.
Thanks

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.

Writing a solution for anyone interested. The issue stems from a dosage error, some variant contains only heterozygous individuals and no homozygous ones needed for QTL calculation.
I used plink2 to calculate genotypes for all variant using the following command:
plink2  --pfile PATH \
     --geno-counts  \
     --out OUT

I then removed variant that contains less than 10 Homozygous invidivuals (HOM_REF_CT < 10) using plink then was able to run tensorQTL interaction mode.
Hope it helps.


text:
TCGA pathology slides magnification
Does anybody know if there is a way to find out what are the magnification levels of the TCGA pathology slides? Are they all 40x? I can't seem to find a label of that anywhere

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.

Can you provide a link to any images? It's likely that the resolution is available in the file metadata.

Thanks. Here is an example of one of them. I can't seem to find the resolution in the metadata

https://portal.gdc.cancer.gov/files/f23db5f3-6b2e-471e-9b18-e477140d9a5e


text:
How do you download metadata for TCGA pathology slides
I downloaded the TCGA pathology slides for the SKCM project and the files are all a name such as "TCGA-3N-A9WB-01Z-00-DX1.A9950ED4-9480-455C-AE0D-8E076D4DA432.svs".
How do you get the clinical metadata for this dataset? When I click "clinical" for the project, the .tsv file has no case IDs that match to my file names.
https://portal.gdc.cancer.gov/exploration?filters=%7B%22content%22%3A%5B%7B%22content%22%3A%7B%22field%22%3A%22cases.case_id%22%2C%22value%22%3A%5B%22set_id%3AYavV3YoBvzAafV-N5zuW%22%5D%7D%2C%22op%22%3A%22IN%22%7D%5D%2C%22op%22%3A%22AND%22%7D

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.

The patient ID is the first 10 letters, e.g. "TCGA-3N-A9WB-01Z-00-DX1.A9950ED4-9480-455C-AE0D-8E076D4DA432.svs" relates to the patient "TCGA-3N-A9WB". If you download the clinical information for SKCM, there is clinical info in the tsv for patient "TCGA-3N-A9WB". Does that make sense?

Thank you! That makes sense!


text:
Ramachandran Plot for specific amino acids in a pdb file
Does anyone know of a Ramachandran plot server (URL please) that allows you to display, for example, just the leucine residues in a protein structure rather than displaying all the amino acids in the pdb file? Thanks!
Paul

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.

https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/ramachandran.html  : Usage: ramachandran [ assign ]  atom-spec
The ramachandran command shows the distribution of peptide φ,ψ angles (a Ramachandran plot) for each specified model that contains protein. It is the command-line implementation of the Ramachandran plot function in the Model Panel. A separate plot will be generated for each protein-containing model with any atom in atom-spec, and the plot for a model will contain all of its peptide φ,ψ angles regardless of which atoms were specified.


text:
How to interpret a GenomeScope plot?
I am working on a new sponge genome and I have produced the genomescope plot from the 21-mer kmer frequencies. I have hard time interpreting the plot. Can someone please help me?
Thank you.

Hi @MudithEkanayake IMO the question has been answered. Please ask via comments if you're not clear about the answer. Do remember to upvote/accept when you're happy with the answer.

Based on the numbers shown above the plot, and the appearance of the curve, this looks like reads from a diploid genome approximately 113Mb in length that was sequenced at around 174X coverage.


text:
MMSeqs taxonomy running for over a day
I've been trying to run mmseqs2 on a few metagenomic assemblies and despite my best efforts in reading the wiki and playing with parameters, the process is taking over a day.
In their paper they claim to be 2 - 18x faster than CAT. I ran CAT on the same datasets and they took at most a few hours. So I can only assume that the long run time for mmseqs is down to user error (i.e. me).
I am currently running the commands (general format):
mmseqs createdb metagenome1_scaffolds.fasta metagenome1_query_DB
mmseqs taxonomy metagenome1_query_DB /home/mmseqs2/NR/NR metagenome1_mmseqs_taxonomy.txt metagenome1_temp --tax-lineage 1 --orf-filter 1 --threads 24

Could someone please let me know where I'm going wrong or what parameters I can tweak to get this running in the expected time frame?

Have you contacted the developers about this? Do any of these issues match the problem you're having?

Hi there! I actually got in running in ~4 hours but then the results are split into 24 semi-corrupted files (I used 24 threads). I have since contacted the developers and will post a solution once I get a response and figure it out :)

After much googling about I realised that I was missing a step! Here are the commands I used to get it running:
mmseqs createdb metagenome1_scaffolds.fasta metagenome1_query_DB
mmseqs taxonomy metagenome1_query_DB /home/mmseqs2/NR/NR metagenome1_mmseqs_taxonomy.txt metagenome1_temp --tax-lineage 1 --orf-filter 1 --threads 24 --split-memory-limit 500G
mmseqs createtsv metagenome1_DB metagenome1_mmseqs_taxonomy.txt metagenome1_mmseqs_taxonomy_final.txt

The split-memory-limit flag sped up the process and the final line put all the split files back together

Thanks @Rainman


text:
Fastest way to convert GTF or GFF file to a BED file?
Given a GTF/GFF3 annotation file, what is the fastest way to convert it into a BED file?

Hi @AlejandroGonzalesIrribarren. Thanks for your response. For future reference please provide some context to this question, e.g. the genome (humans), computer resources and what have your tried? Again thanks for posting back. Conda is cool for getting stuff working quickly, but not necessarily performant .... Easybuild is much better in this particularly question, i.e. fastest.

Yes, please at least provide example input and expected output so the question can be clear.

gxf2bed
The fastest G{F,T}F-to-BED converter around the block!
translates
chr27 gxf2bed gene 17266470 17285418 . + . gene_id "ENSG00000151743";
chr27 gxf2bed transcript 17266470 17281218 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8";
chr27 gxf2bed exon 17266470 17266572 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8"; exon_number "1"; exon_id "ENST00000541931.8.1";

into
chr27 17266469 17281218 ENST00000541931.8 1000 + 17266469 17281218 0,0,200 2 103,74, 0,14675,

before your eyes blink!
Converts

Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 3.34 seconds.
Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 2.11 seconds.
Canis lupus familiaris ROS_Cfam_1.0 Ensembl 110 (55,335 transcripts) in 0.99 seconds.
Gallus galus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.14 seconds.

Usage
Usage: gxf2bed[EXE] --input/-i <GTF/GFF> --output/-o <BED>

Arguments:
    --input/-i <GTF/GFF>: a .gtf/.gff file
    --output/-o <BED>: path to output .bed file

Options:
    --help: print help
    --version: print version
    --threads/-t: number of threads (default: max ncpus)

Installation
to install gxf2bed on your system follow this steps:

get rust: curl https://sh.rustup.rs -sSf | sh on unix, or go here for other options
run cargo install gxf2bed (make sure ~/.cargo/bin is in your $PATH before running it)
use gxf2bed with the required arguments
enjoy!

Build
to build gxf2bed from this repo, do:

get rust (as described above)
run git clone https://github.com/alejandrogzi/gxf2bed.git && cd gxf2bed
run cargo run --release -- -i <GTF/GFF> -o <BED>

Container image
to build the development container image:

run git clone https://github.com/alejandrogzi/gxf2bed.git && cd gxf2bed
initialize docker with start docker or systemctl start docker
build the image docker image build --tag gxf2bed .
run docker run --rm -v "[dir_where_your_gtf_is]:/dir" gxf2bed -i /dir/<GTF/GFF> -o /dir/<BED>

Conda
to use gxf2bed through Conda just:

conda install gxf2bed -c bioconda or conda create -n gxf2bed -c bioconda gxf2bed

Benchmark + FAQ

If you google "gtf to bed" or "gff to bed" you'll find some posts about recommended tools. Part of these tools are already deprecated or too difficult to run (even in Linux), due to poor maintenance. This project was conceived to provide an easy way to convert GTF/GFF files to BED structures and finish a set of high-performance converters between gene model formats in Rust (bed2gtf, bed2gff and now gxf2bed).
To test the efficiency of gxf2bed, I took 4 GTF/GFF-to-BED converters and ran them with the current GRCh38 GENCODE annotation (v.44), that has ~250,000 transcripts [the values showed here are the mean of 5 consecutive runs]. Briefly, the competitors are:

UCSC's utils: gtfToGenePred | genePredToBed & gffToGenePred | genePredToBed (taken from here)
Signal & Brown: gtf2bed.py & gff32gtf.py && gtf2bed.py (taken from here)
ea-utils: gtf2bed.pl & gff2gtf.pl | gtf2bed.pl (taken from here)
BEDOPS: gtf2bed &  gtf2bed (taken from here)

Besides the easy way to make gxf2bed run and the GFF-GTF channel that allows centralize both in a single converter, this tool showed a significant difference in computation time against the other tools in each one of the two formats. Even using a combination of different number of threads, times were practically the same 3.5s +/- 0.4s and the differences were maintained. Is important to note that GFF3 files took a lot more time than its GTF counterpart. This could be due to the need of chaining two different programs 1) GFF-to-GTF and 2) GTF-to-BED.
Taken together, gxf2bed offers an easier and faster way to convert GTF/GFF files to BED files. This tool could save you at least x2-3 times the time you used to spend using other tools to convert GTF files and x5 times if you want to convert GFF3 files!

Disclaimer: I am the author of the gxf2bed tool.

How do you deal with regular GTFs that have a variable number of key-value pairs in the ninth column?

It does not matter how many number of key-value pairs are in the ninth column. The only thing this tool assumes is that at least "gene_id" and "transcript_id" should be present in each line. This is nothing crazy to ask since is very common practice.

So what options did you use to make the conversion you show here? Why were the gene and exon lines ignored and only the transcript kept?

Exactly what @terdon says - there doesn't seem to be any way to pick the features to extract. Is your tool faster than a simple awk in that case? And are you sure you're comparing apples to apples when your tool makes these assumptions?

Exon lines are not ignored. The tool basically parses a gtf and extracts only the needed features to build a BED12 (which is basically a transcript-based annotation file). The results of gxf2bed are exactly the same obtained by UCSC's binaries and the other tools as well. gxf2bed even goes a step further sorting the BED12 by chr and pos, which does not seem to be accurate for some of the tools evaluated.
If you are confident that an awk command could build a correct BED12, wonderful! This tool provides another opportunity to perform this conversion in less time that the available tools now.

If you are skeptical about how the results of gxf2bed, you can try it, though. Just run the command as explained in this post and compare the results yourself. To elaborate a little bit more on the last lines of the past comment: re-writting this functionality in Rust now allows 1) a more straightforward installation, 2) faster times, 3) easy pipeline implementation, 4) better performance. All this is open source, so feel free to review the code if you still do not feel that comfortable! :)

"a BED12 (which is basically a transcript-based annotation file)" Please show your source for this information. BED files are lists of contigs+positions with optional annotations, they could be genes, transcripts or even random nucleotides. You are assuming the transcript part.

gtf2bed, for example, accepts an attribute_key and is thus flexible. gtf2genePred does assume gene as the feature and it does provide some flexibility as to the identifier within that constraint, but it calls itself genePred, not featurePred. Sure, gtf2bed.pl seems to default to transcript_id but that's the only tool you can compare to. I appreciate your rust-based contributions to the field but please be more mindful with your claims.


text:
How to interpret high support values for terminal nodes and low for more ancient nodes?
I made a phylogenetic tree for about 300 species using a region of the COX1 gene. The tree looks visually meaningful: species of the same genus are all clustered together.
Then I performed a 350-replicate bootstrapping and I analyzed the support values. Surprisingly I found very low values for nodes that are more close to the root and high support values for nodes that are more close to the leaves.
Almost all nodes that split into 2 leaves have Bootstrap >95% but the ones close to the root have Bootstrap<20%.
Does that mean I chose a highly variable DNA region?
How would you re-build the tree to make it more consistent on more ancient nodes? I was thinking of using a concatenation of sequences with both high and low variability regions, but i am not sure it would work.

What are the families and the Order under investigation? Is this a nucleotide or protein tree?  If you are considering the 'control region of mtDNA it suggests this is a nucleotide tree.

This is the expected result. We can be much more certain about the terminal nodes of the tree than we can be about the tree's deep structure, trivially (but not exclusively) because of uncertainty propagation from the tips down into the internal nodes (a deep node is only as good as all of the nodes that lead to it). You can see a similar phenomenon in the various phylogenies in this classic paper.
It is also possible that your gene is not variable enough, not just that it is too variable. If there is a subset of highly informative positions in your gene, bootstrapping samples that miss those positions would be expected to yield noisy results.
In terms of how to resolve these, the traditional answer is "more data". This data can come in either the form of more species/tips or more sequences per species (ideally both). It is probably unsurprising that you aren't able to get very confident with a single gene. Even a rather older approach uses 31 universally conserved genes. As is traditional in statistics, we get more certain about our estimates as we use more data to arrive at them. In these methods, the genes (proteins) are indeed concatenated.
Note also that single gene trees can be positively misleading about the true phylogenetic relationship, due to incomplete lineage sorting or other such phenomena. In this case, you could have excellent bootstrap support for a wrong phylogeny.
In another paper, you will find the following statement:

In molecular phylogenetic analysis, it is often the case that some relationships are robust, and relatively "easy" to recover while others are difficult to resolve, leading to phylogenetic hypotheses that consist of a patchwork of well and poorly supported nodes. When difficult nodes are encountered, the next logical step is to add taxa and/or data under the reasonable assumption that additional taxa or characters might enable resolution and/or provide support for poorly supported nodes...

The amount of data required for resolution of difficult phylogenetic problems associated with short internodes, especially those deep in a tree can represent a particularly difficult challenge [5–7] that often requires massive amounts of sequence data to resolve.

Deep nodes are just harder than shallow nodes! The possible exception here is the very oldest splits, where enormous amounts of evolutionary time might separate outgroups quite clearly from a relatively messy subtree (see the Fishbein paper for more examples of this).


text:
BioPython bootstrap is not reliable?
Here i will show you a minimal working example of code and as you can see the support values for the tree is always 100.
I am using synthetic sequences of 100bp for 6 elements. The sequences have been generated at random choosing from ATCG for each position with equal probabilities.
The input fasta file is
>S1
CTCAGGCACTAGGAGCTTCTCCAGGGCAAAGTTGTCTACGAATATGCCGACTCAGAAGGTTATCAATACGGTTTACTTATCTGCACGCCATTTTCCTATG
>S2
TAACAAATGTTTCTCGCTCGAACCCCGTGGTGCGAGGGGTCACGATAAAAGGTCCCTCCTTCCACGGCATAATTGCTCCCTTTCTTTTCCGTGGGCAGGA
>S3
TTCCTAAAGGACGACTGGAAGCCGGGTACACGTCAACGGAGCTTCTAGCCCGAGGTCTATAGACCGGTTAATGAACTAGACTTAATGTGGAGCTCGTAGA
>S4
ACTTCGTGCCCAACCACTCTATGAGCAGAGTGCTCGAATGAACCTCAAAAGGATATCCGCATTTACTCTATATAACAAACGGCCGTCCCCCCATCTGTCA
>S5
TGACACTGGGTCATTTTACCCGCACTGATCGCTGGGCAGGTCGGCAATTTCGTCAGAATGCCGTGGCCGCACTGAAAATATCATTACCGGACTGAGTATG
>S6
ATAAGCCGAGGGGTAGCCCTCTATTTCGCACCGTATAGACGAAGTGATAAACTTTCTAACCACGTGTCGCCTATCTCTACCTAGCCACATTTGAGGTGCG

The following python code makes the multiple alignment using MUSCLE algorithm
### alignment
from Bio.Align.Applications import MuscleCommandline
from Bio import AlignIO, SeqIO
import os

muscle_exe = "C:/Scripts/muscle3.exe"
path = os.getcwd()
with open('tmp_in.fa', 'w') as fp:
    fp.write(fasta)

## replacing \\ with /
inp_fp = path.replace('\\', '/') + '/tmp_in.fa'
out_fp = path.replace('\\', '/') + '/tmp_out.afa'

muscle_cline = MuscleCommandline(muscle_exe, input=inp_fp, out=out_fp)

The following code block takes the previous alignment out_fp and creates a tree using UPGMA algorithm.
The number of bootstrap trees is set to 50 and the consensus is majority based.
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo.Consensus import *

# read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')

calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)

## get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)

print(consensus_tree)

The final consensus tree looks like this
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.350925925925926, name='S5')
        Clade(branch_length=0.05833333333333341, confidence=100.0)
            Clade(branch_length=0.05879629629629617, confidence=100.0)
                Clade(branch_length=0.26481481481481484, name='S3')
                Clade(branch_length=0.020370370370370372, confidence=100.0)
                    Clade(branch_length=0.24444444444444446, name='S1')
                    Clade(branch_length=0.24444444444444446, name='S4')
            Clade(branch_length=0.031018518518518404, confidence=100.0)
                Clade(branch_length=0.2925925925925926, name='S6')
                Clade(branch_length=0.2925925925925926, name='S2')

My question is: why does the tree have all support values equal to 100? I am sure there must be an error either on my side or in the package module itself.

I think this is a bug.
It seems to work if you do this, creating an equivalent Alignment object instead of a MultipleSeqAlignment to give the bootstrap step:
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)

bootstrap_consensus calls bootstrap_trees, which makes however many randomly shuffled alignments you asked for.  But if the input is a MultipleSeqAlignment object like you're using it builds a tree over and over with that same input to the function as the tree function input, so it always gets the same tree over and over.  The other branch of bootstrap_trees, which runs if you have an Aignment object, builds trees based on the randomly-selected columns of the input, like you'd expect.  Pretty sure that's just a variable naming screw-up one of us should report to Biopython.
Self-contained version of your tree part:
#!/usr/bin/env python
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio import AlignIO
from Bio.Phylo.Consensus import *

out_fp = "tmp_out.fa"

# read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')

calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)

## get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)

print(consensus_tree)

# Give bootstrap_consensus an Alignment object instead of MultipleSeqAlignment
# object
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)
print(consensus_tree)

Example outputs for the two different input alignment objects:
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.012239583333333349, confidence=100.0)
            Clade(branch_length=0.31666666666666665, name='S3')
            Clade(branch_length=0.033333333333333284, confidence=100.0)
                Clade(branch_length=0.2833333333333333, name='S2')
                Clade(branch_length=0.2833333333333333, name='S4')
        Clade(branch_length=0.033072916666666716, confidence=100.0)
            Clade(branch_length=0.29583333333333334, name='S5')
            Clade(branch_length=0.03750000000000007, confidence=100.0)
                Clade(branch_length=0.2583333333333333, name='S1')
                Clade(branch_length=0.2583333333333333, name='S6')
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.2583333333333333, name='S2')
        Clade(branch_length=0.03405145202020203, confidence=66.0)
            Clade(branch_length=0.275, name='S5')
            Clade(branch_length=0.049892241379310344, confidence=58.0)
                Clade(branch_length=0.21250000000000002, name='S1')
                Clade(branch_length=0.21250000000000002, name='S6')
        Clade(branch_length=0.039416666666666676, confidence=50.0)
            Clade(branch_length=0.2583333333333333, name='S4')
            Clade(branch_length=0.321875, name='S3')

Thank you very much i think this was the issue. I will report it to their github if you haven't already raised an issue.

I haven't yet, so, go for it!  Should probably be a simple fix for them.

Initially, 100 looks like a good value.  However, if you lower to bootstrap value, let's say maybe to 50 as you did, you will probably get an erroneous support value.  Bootstrapping (you probably already know this) is how many times a branch in a phylogenetic tree is observed.  A good bootstrapping value would be anywhere between 80 and 100.  50 is not a good bootstrapping value.  This means the branches are only being observed 50 times as opposed to 85, 90, etc...   If I were you, try changing the bootstrapping value of 50 to at least 80.  See what that gets you.  You may be satisfied with it, but you may want to increase it depending on your output.  By doing this, you will see a more realistic and phylogenetically accurate output of your tree.  Most bioinformaticians use 100 as a bootstrap value.  The only issue, I wouldn't really call it an issue, but it will take longer for an output to occur due to computing time and computing memory available.  BUT you will still get an output, you just may need to wait a bit.   My bioinformatics background is rooted in phylogenetic and horizontal gene transfer research involving phages and bacteria, so I have seen this issue a lot and have experienced myself a lot!  I hope this information provides some insight to your project and bioinformatics journey! :)

I don't think this is the case. This is a Minimal working example to show the problem. The real data is hundreds of species with longer sequences. When I run multiple times manually the clades in the tree change, but when I use biopython consensus tree it does not happen, so it must definitely be a bug or an error in my code.

When you mean you ran the samples "multiple times manually" as opposed to "when I use the biopython consensus tree"?  This sentence is confusing and misleading potentially.  "Consensus trees are used to summarize a set of trees defined on the same set of taxa. A consensus algorithm takes the trees as input, so that the method of producing the input trees is not part of the consensus algorithm" as cited on:  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2909780/#:~:text=Consensus%20trees%20are%20used%20to,part%20of%20the%20consensus%20algorithm.

A consensus tree uses the consensus algorithm to generate a tree.   That makes me wonder if the algorithm you are using is appropriate for the amount of samples you have and the variation of species from different axons that you are using.  What taxons are you using?  Are you looking at family, class, order, genus, etc...?  Family would defiantly be a broader range of values than genus.  If you can maybe narrow down to a more exclusive taxon like order or genus, maybe your values may be more consistent?

In the example i gave (that is completely sufficient to understand the problem), there is no mention of taxa. There is no need to invoke more complex processes when it is evident that when you run a consensus tree either on a very random sequence set (MWE) or in a very big dataset (my actual case) you do not expect to see 100% on all the nodes. It's like saying that you made the best possible alignment, also all the assumptions of your models are met and there is no measure error in the sequences and so on...


text:
Is it possible to, all in memory without writing a file, given a list of AlignedSegments and a bam header, build a bam, index and get coverage?
My current job is rate-limited chiefly by the IO overhead of writing reads to their corresponding bam files (they are split up by both cell barcode and contig). The only reason I need to write these mini-split bams is so that I can then use the pysam coverage function to assess read depth at certain locations. I already have all of these reads in memory as a list of AlignedSegments. Is it necessary to write them to a .bam file in order to be able to index and use the coverage function? I feel like there should be a way to have the AlignedSegments be treated as a bam object in memory, possibly even indexed in memory, and then treated as if they were just streamed from a bamfile with an index in order to be able to use the coverage method. Thanks in advance, if this is possible it would increase my processing speed x10. Or, given a list of AlignedSegments, is there a way to get coverage at specific positions across all those AlignedSegments in my list (that would get me what I need).

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.

What operating system are you using?

My operating system is linux

Assuming you are on an operating system that supports them, like Linux, you can try using named pipes in a ramdisk or tmpfs. Not the most elegant approach, but I don't know if what you are asking for is possible (it may very well be, I just don't know) and this can be a decent workaround.
First, check that you do have a mounted tmpfs. On my system, I have a few:
$ mount | awk '$1=="tmpfs"'
tmpfs on /dev/shm type tmpfs (rw,nosuid,nodev,inode64)
tmpfs on /tmp type tmpfs (rw,nosuid,nodev,size=16278892k,nr_inodes=1048576,inode64)
tmpfs on /run/user/1000 type tmpfs (rw,nosuid,nodev,relatime,size=3255776k,nr_inodes=813944,mode=700,uid=1000,gid=1000,inode64)

These filesystems are actually stored in RAM and should be significantly faster. For instance, on my machine, which has an SSD as the main storage, I still get significant speedup when using /dev/shm:
## SSD
$ time seq 1000000000 > file
real    0m13.687s
user    0m5.623s
sys     0m7.516s

$ time wc -l < file
1000000000

real    0m4.972s
user    0m0.101s
sys     0m2.248s

## tmpfs, RAM
$ time seq 1000000000 > /dev/shm/file

real    0m6.755s
user    0m4.982s
sys     0m1.763s

$ time wc -l < /dev/shm/file
1000000000

real    0m1.024s
user    0m0.120s
sys     0m0.902s

As you can see, using /dev/shm is significantly faster. And using named pipes can give you a little boost too:
$ time ( seq 1000000000 > /dev/shm/file;  wc -l < /dev/shm/file)
1000000000

real    0m8.341s
user    0m5.277s
sys     0m3.044s

$ rm /dev/shm/file && mkfifo /dev/shm/file
$ time ( seq 1000000000 > /dev/shm/file &  wc -l < /dev/shm/file)
1000000000

real    0m6.564s
user    0m5.426s
sys     0m3.965s

So, this is essentially a way of doing everything in RAM without needing to code it as doing everything in RAM. You can create these files (the FIFOs, the named pipes) and treat them as files, but since you're creating them in a tmpfs filesystem, they're actually in RAM and by using FIFOs you not only get a further speed jump but also minimise the actual data that need to be stored.

This seems like it gets at what I'm trying to accomplish. I will try it! Thank you.