text:
Ubiquitous regulation of highly specific marker genes
I am fairly new to scRNAseq analysis and keep running into the same problem in the two datasets I am currently working with. We work with kidney inflammation and have two new mouse models, which we are now analysing in scRNAseq. I followed the usual Seurat workflow and ended up getting some nice UMAPs. With FindAllMarkers I was able to identify the clusters pretty well.
But once I do differential expression analysis (either with FindMarkers or pseudobulking) I got weird results.
Background There is a highly specialized epithelial celltype in the kidney (the podocyte) that specifically expresses a host of markers (like e.g. Nphs1 and Nphs2). They get downregulated when the kidney gets damaged (as it does in our disease models).
Problem

I do not expect to see those markers to be present or regulated in any other celltype. But when I do the DE, those podocyte-specific marker genes pop up in almost all other clusters and are described to be highly regulated (downregulated, as I would expect in the pdoocytes).

So in my eyes, this points to "doublets of podocytes" and other celltypes being involved here. Maybe some damaged podocytes got into the same droplet as other celltypes and thus distort the analysis.

I then used somewhat stringent filtering via scDblFinder* to remove a lot of doublets (~ 10 % of my dataset).

However, I still get the same problem.

Then I analysed each batch (we have 3 in total) individually, but in each one this problem comes up.

Overall, I see 6.4% of doublets (that is 8,710 doublets in a total of 135,844 cells). I have two podocyte clusters with 1.7% doublets (that is 462 doublets in a total of 26,693 podocytes).
Some help / advice would be highly appreciated as I am becoming a little frustrated with this issue.
Thanks a lot!
*, Bioconductor package

An explanation in 3 parts
Final Edit:  It turns out there was contamination in the data which you can read about in the comments, without being able to see it myself I could not know this.  If you see this in an experiment and you're clustering from totally raw unfiltered data this maybe the case for you.  Nevertheless as you subcluster within cell types of particular interest to get more refined marker genes the advice I provide below may be of use.
EDIT: The note above about looking at the violin plots may be useful, but again what your downstream analysis or experiment should impact your decision making.
TLDR:  The presence of marker genes in other cell types is not surprising nor necessarily a problem.  A marker gene that seurat finds are highly expressed in a cell type relative to the average in all cells not a gene that is uniquely expressed, if you are looking for uniquely expressed genes you will probably not find them.
1) A note on seurat
So most genes are expressed in most cells there is nothing wrong with your data or your analysis.  Gene expression is very messy and the term "marker gene" is a bit of a misnomer.  Note that other analysis pipelines like Monocle have different methods for determining what constitutes a marker gene.
Please consider this meta-analysis of gabanergic cerebral cortex neurons (i.e. inhibitory cortical neurons, interneurons etc etc).
https://ieeexplore.ieee.org/abstract/document/9669888
The fact is that most cell populations express a gradient of genes and that each gene likely has multiple roles in each cell.  If your goal is to simply publish on expressed marker genes taking the top highest uniquely expressed genes which achieve significance in each cluster will be sufficient, just be sure to report ALL your settings and provide the meta data when you publish.  If you are going on to do MERFISH or stain with the results you may want to be more selective.
FindaAllMarkers() uses a simple DEG analysis where it takes the average expression of each feature (the genes you selected, which I hope you did because it's not wroth looking at every gene for clustering just the top most variable) in the cluster against the average expression for these genes in all other clusters.  Note that this means you are doing A TON of statistical tests so they use the bonferroni correction which would normally be considered to conservative but in this case is probably still too lenient.  You should only be reporting on the very top 5 maybe 10 expressed genes. If your lucky one of these will be noticeably more prominent than the others.
2) Alternatives to seurat
Monocle package has a different definition which arguably is better: http://cole-trapnell-lab.github.io/monocle-release/docs/#classifying-and-counting-cells though this package is less widely used.  The monocle score accounts for how widely expressed across all clusters a marker is and the more unique to a cluster a gene is the greater its score in addition to its actual expression strength.
There is also a package called pandora, but I know next to nothing about it.  Read up on them and decide which has features you like best.
Lastly I think part of your problem is just wanting biology to be a cleaner science than it is.  The reality is that things we call cell types are simply transient states a population of cells passes through during development and various environmental perturbations, the reality is each cell is a dynamic system and you are sampling a portion of these graded changes.  Some states are more stable than others and we have taken to calling these cell types but remember you can always subcluster within a cluster and pull out more clusters which can themselves be subclustered and so on and so forth.  It calls into question the definition of a cell type even more and to be honest at this point finding cell types is somewhat of an art, the science is yet to be settled.
I encourage you to do more reading on gene expression within single cells and the definition of a cell type (or if such a thing even exists).  Here is some reading to get you started.
Cembrowski, Mark S., and Vilas Menon. "Continuous variation within cell types of the nervous system." Trends in Neurosciences 41, no. 6 (2018): 337-348.
https://www.sciencedirect.com/science/article/pii/S0166223618300596
Muñoz-Manchado, A. B., Gonzales, C. B., Zeisel, A., Munguba, H., Bekkouche, B., Skene, N. G., ... & Hjerling-Leffler, J. (2018). Diversity of interneurons in the dorsal striatum revealed by single-cell RNA sequencing and PatchSeq. Cell reports, 24(8), 2179-2190.
https://www.sciencedirect.com/science/article/pii/S2211124718311562
If you really want to lean into this you should start digging into the concept of meta-cells, there is a package called Metacell2 which may be of use, be sure to read the documentation and methods section for Metacell (the first version) as Metacell2 won't make much sense without it:
Harris, Benjamin D., Megan Crow, Stephan Fischer, and Jesse Gillis. "Single-cell co-expression analysis reveals that transcriptional modules are shared across cell types in the brain." Cell Systems 12, no. 7 (2021): 748-756.
-this is not metacell but is a meta analysis of cell markers using metacells with some pretty surprising and counter intuitive results about marker genes
Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., ... & Tanay, A. (2019). MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions. Genome biology, 20(1), 1-19.
Ensure you read all documentation carefully and understand the stats they are using.  The methods use are a specific definition.  Also a real marker gene will be robust to multiple methods of identifying them, if they disappear based on a change in normalization or in other packages you may have your answer there.  Bear in mind if you publish your results should be reproducible and people will need to know exactly how you got your results so share your code and metadata please.
3) What do you want to get from this analysis?
As I mentioned earlier if your goal is to get published, simply finding some markers in seurat may suffice.  If your goal is to do further experiments based on these results consider your hypothesis carefully.
If you want more certainty pull some similar kidney scRNA seq data from GEO and compare your results to their or integrate your data with theirs seurat has a tutorial you can follow along with to do this https://satijalab.org/seurat/articles/integration_introduction.html#performing-integration-on-datasets-normalized-with-sctransform-1 Note you can also use this pipeline to do batch correction, that may help you.
The sct transform give pretty good results.
Sorry nature isn't like a biology textbook =(,  I'm not saying that to be facetious, I just want to remind you that living things are very complex.  Lean into that complexity and make some interesting discoveries or just accept it and dig out some results you can publish.  Talk with your supervisor and explain the issue.  Hope this helps!
4) I lied sorry this is a response in 4 parts, make sure you report exactly what you did and share your metadata and pipeline
Did I mention your results should be reproducible?  How can one do this?  Here are some tips.  Make sure your keeping your analysis on github, so you can keep track of all the changes it goes through and revert to previous versions if you change things.  Keep all your code, its worth having a file that is full of old code you threw away cause it was messy or you didn't end up using it. Not broken code but stuff you don't think is useful.  Report your final results, share your data and make sure a monkey could follow you instructions to reproduce your results as well as being able to double check they get the same thing.  Use a conda environment to ensure software updates will not destroy your code for future lab members or people who read your future paper.  Save all file paths as variables at the top of your code in your scripts so people can easily change the file structure without having to go through your code and struggle finding which read.csv() functions they need to edit.

Thanks so much for your detailed answer Angus, this was really very helpful for me as a beginner. As mentioned in my answer below, we discussed the dataset with more experienced bioinformaticians and it seems like there was a big amount of ambient RNA contamination in the dataset. Correction for this with the DecontX package cleaned things up a good bit (even though of course there is still some celltype-specific transcripts floating around in the other celltypes which is expected to a certain degree, as you mentioned).

Glad you caught that confound.  I work in brain tissue, which is much messier often people are doing lots of sub clustering, i.e. they do the initial round of clustering but then select only the neurons and redo the process.  Then they might pull out only the excitatory or inhibitory cells (the two major types of neurons) and subcluster within just those cells again.  As you imagine it gets messy.  I can see how if at the initial round of clustering you are seeing this extreme level of non-specificity that would suggest contamination.  Glad you caught it.

Hi @Smeerlap could you please very kindly, upvote and/or mark this answer as "accepted". It is a very cool answer as you've said.  Seurat vs. Monocle was very helpful for me and first time I've seen a comparison on this site.

The violin plots show only little expression in non-podocyte clusters, but the podocyte marker transcripts are only the most obvious manifestation of this problem. Also, small but significant amounts of endothelial and other celltype transcripts show up in most of the other clusters. Identification of the celltype works fine but these contaminating counts throw off most of the downstream analyses.
By now, we did manage to find a solution for the problem, though: It turns out the problem was due to high amounts of ambient RNA contamination. I tried three different packages to remove the contaminating counts (DIEM, DecontX and SoupX) and found that DecontX (and SoupX) worked best for my dataset.

You should make violin plots of the genes of interest (x-axis is your different clusters (i.e. cell-types) and y-axis is expression). Maybe they are differentially expressed but at the same time not expressed very highly.
Also why does their expression matter in other cell types? Do you have strong evidence to expect that these genes would only be expressed in podocytes?

Hi! Thanks for your help. It seems like the problem was due to ambient RNA contamination (see my post below).


text:
When does a 2D histogram / heatmap take a circular shape?

The above is a 2D histogram of y and z internal coordinates describing the positions of C-alpha atoms involved in hydrogen bonds between antiparallel Beta-strands.
There is no color coding. It is given as it is.
What can we assume about the data points if a histogram is circular?
When does a histogram take a circular shape?
[Note: "2D histogram" is a synonym for what bioinformaticians typically call heat maps, see here]

source: Stack Exchange

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

sort: votes


text:
What does it mean by a Local Coordinate in the case of protein chains?

High-performance transformation of protein structure representation from internal to Cartesian coordinates

In Figure 3a, each residue has its own local coordinates, which are illustrated shown using different colors.

Here are two more papers:

Generative De Novo Protein Design with Global Context

Generating 3D Molecules for Target Protein Binding

By global coordinate, I understand the coordinate system in which we view the entire protein chain.
What does it generally (not only in this particular paper) mean by a Local Coordinate in the case of protein chains?
Why was it necessary to introduce this terminology when they could have just used the term "vector"?
Can you supply me with an example with an explanatory diagram?

For proteins and molecules in general, and within the meaning of this paper, "internal coordinates" are any representation that are not dependent on the molecule rotation or translation in a Cartesian space.
Calling them "internal coordinates" may be slightly confusing, and a more correct term could be "internal representation"
That means, I guess, that internal coordinates could be, as you mention, a list of vector in a Cartesian space, where the Cartesian space "orientation" would be decided by, say, principle component analysis (or any other arbitrary rule depending on the molecule). This however defeats the purpose of most internal representations: dimensionality reduction or faster computations.
Other (more interesting) ways of representing a molecule in internal coordinates would be:

As (Bayati, Lardan & Bardhan 2020) focus on: a set of bond distance, angle, and dihedrals (see Ramachandran plots and rotamers). This is especially useful in the case of proteins. This representation is for example used by Rosetta. Angles and bond distance can be considered as constant (depending on the atomic type), and so one can infer the Cartesian coordinates of multiple atoms with just one dihedral angle, hence providing a huge dimensionality reduction. In this case there is no need for vector. Also, in the case of internal energy computation there is always one energy term depending on dihedrals, so that is already done. Limit of such methods is that in molecular constructs there are energy terms due to long-range interactions, i.e van der waals and electrostatics. There is hence a need in such software to always switch between a Dihedral representation and a Cartesian representation. Another problem is that proteins backbone are then represented as a dihedral "chain", and so the Cartesian coordinates of one atom becomes tightly linked to the previous one, which are linked to the previous one, etc. This is the problem this paper is trying to deal with.

As (Tan et al. 2023) use, a graph representation. Graphein is an interesting python library for converting and analyzing proteins and molecules from cartesian coordinates to graph representation. In this case instead of having a vector, each atoms position is defined as its distance with surrounding atoms. This kind of representation has received more and more attention from the protein structure prediction & design community as it is particularly suited for neural networks. For example I would say the protein representation in AlphaFold is a kind of graph representation (using "triangles" made by the distance between three atoms). In comparison to dihedral representation, graph representation packs less "assumption" about the system and are even worst than Cartesian representation in terms of dimension (if you do a "global" representation of N atoms then you end up with N! values vs only 3*N with Cartesian), but at least things like non-covalent interactions are implicitly represented.

I believe this answer gives enough material to look at. If this still appears unclear, skimming through Rosetta tutorials / papers figures / webinars should provide ample explanations, same for AlphaFold.

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


text:
Error in dimnamesGets(x, value) when trying to read data using Seurat package
This question has also been asked on Biostars
I am trying to create a Seurat object using the package "Seurat". When I am reading my raw files using the function Read10X, I am getting this error
Error in dimnamesGets(x, value) : length of Dimnames[[1]] (45948) is not equal to Dim[1] (45947).

I am unsure how to resolve this issue.
R sessionInfo:
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sp_1.5-0           SeuratObject_4.1.1 Seurat_4.1.1       Matrix_1.5-1

loaded via a namespace (and not attached):
  [1] nlme_3.1-153          spatstat.sparse_2.1-1 matrixStats_0.62.0    RcppAnnoy_0.0.19
  [5] RColorBrewer_1.1-3    httr_1.4.4            sctransform_0.3.4     tools_4.1.2
  [9] utf8_1.2.2            R6_2.5.1              irlba_2.3.5           rpart_4.1-15
 [13] KernSmooth_2.23-20    uwot_0.1.14           mgcv_1.8-38           rgeos_0.5-9
 [17] DBI_1.1.3             lazyeval_0.2.2        colorspace_2.0-3      tidyselect_1.1.2
 [21] gridExtra_2.3         compiler_4.1.2        progressr_0.11.0      cli_3.4.0
 [25] plotly_4.10.0         scales_1.2.1          lmtest_0.9-40         spatstat.data_2.2-0
 [29] ggridges_0.5.3        pbapply_1.5-0         goftest_1.2-3         stringr_1.4.1
 [33] digest_0.6.29         spatstat.utils_2.3-1  pkgconfig_2.0.3       htmltools_0.5.3
 [37] parallelly_1.32.1     fastmap_1.1.0         htmlwidgets_1.5.4     rlang_1.0.5
 [41] rstudioapi_0.14       shiny_1.7.2           generics_0.1.3        zoo_1.8-10
 [45] jsonlite_1.8.0        spatstat.random_2.2-0 ica_1.0-3             dplyr_1.0.10
 [49] magrittr_2.0.3        patchwork_1.1.2       Rcpp_1.0.9            munsell_0.5.0
 [53] fansi_1.0.3           abind_1.4-5           reticulate_1.26       lifecycle_1.0.2
 [57] stringi_1.7.6         MASS_7.3-54           Rtsne_0.16            plyr_1.8.7
 [61] grid_4.1.2            parallel_4.1.2        listenv_0.8.0         promises_1.2.0.1
 [65] ggrepel_0.9.1         deldir_1.0-6          miniUI_0.1.1.1        lattice_0.20-45
 [69] cowplot_1.1.1         splines_4.1.2         tensor_1.5            pillar_1.8.1
 [73] igraph_1.3.4          spatstat.geom_2.4-0   future.apply_1.9.1    reshape2_1.4.4
 [77] codetools_0.2-18      leiden_0.4.3          glue_1.6.2            data.table_1.14.2
 [81] png_0.1-7             vctrs_0.4.1           httpuv_1.6.6          polyclip_1.10-0
 [85] gtable_0.3.1          RANN_2.6.1            purrr_0.3.4           spatstat.core_2.4-4
 [89] tidyr_1.2.1           scattermore_0.8       future_1.28.0         assertthat_0.2.1
 [93] ggplot2_3.3.6         mime_0.12             xtable_1.8-4          later_1.3.0
 [97] survival_3.2-13       viridisLite_0.4.1     tibble_3.1.8          cluster_2.1.2
[101] globals_0.16.1        fitdistrplus_1.1-8    ellipsis_0.3.2        ROCR_1.0-11
```

Check the length of your files. The error clearly says that there is a difference in the length of your files. Your file length should be equal.
I think the problem may be the presence of a header in one of your files. Removing the header may solve the problem.

There is a mismatch of 1 line between the several files, suggesting that probably either the barcodes of features file has a header line which is not standard. Remove any header lines from these files.


text:
What is the difference between profile and pattern in Prosite?
I found a definition:
"A pattern is thus a qualitative description of a motif in terms of amino acid sequence. The concept of a profile extends this concept, allowing a quantitative description of a motif, by assigning probabilities to the occurrence of a particular amino acid at each position of a motif."
This is example of profile in prosite
Hits for all PROSITE (release 2022_04) motifs on sequence P98073 [UniProtKB/Swiss-Prot]:
Example of profile
PS50024   SEA   SEA domain profile :
54 - 169:   score = 32.979
QSHEARATFKITsGVTYNPNLQDKLSVDFKVLAFDLQQMIDEIFLSSNLKNEYKNSRVLQ
FEN--GSIIVVFDLFFAQWVSD--ENVKEELIQGLEANKssQLVTFHIDLNSVDILDKLT

This is example of pattern in Prosite
PS01209   LDLRA_1   LDL-receptor class A (LDLRA) domain signature :
197 - 221:     [confidence level: (0)]  CIkadlf.CDgevNCpdgsDEDnkm..C

What is the difference in that case? So in profile I have probabilities (in definition I have "profile extends this concept, allowing a quantitative description of a motif, by assigning probabilities to the occurrence of a particular amino acid". I don't see any probabilities in that sequence.
What does it mean "-" in profiles? Maybe any aminoacid or gap?
Why I have sometimes lowercase letter and uppercase letter in profile and pattern? I don't get it. It's a probability?
What is more important profile or pattern?

Hi @MrMorgotheron and welcome. From my point of view its a good question, but a few too many questions, albeit I could be wrong. The "-" definitely means an indel, so you can remove that question to help focus the question more.

I think this was already answered on a different SE site: https://biology.stackexchange.com/questions/46751/difference-between-motifs-domains-patterns-signatures-and-profiles

This publication explains it nicely: https://academic.oup.com/bib/article/3/3/265/239712?login=false


text:
Benchmarking for variant identification using RNA-seq data
I am in need to benchmark the variant identification pipeline which uses RNA seq data alone without any matched-normal. I would like to know the reference dataset (and the pipeline on which the dataset was run) that can I use to compare the pipeline results on my dataset. The reference dataset which I use should be comparable regardless of the dataset specific to certain disease. Thanks all in advance for your responses.

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.

Well, what is your dataset? What species? What kind of sample? You mention the absence of a matched normal, does that mean your sample is from a tumor? We need to understand what you have in order to help you validate.

The samples come from human related to diabetic condition


text:
Allele effects by parental-origin estimated from genotyping array data
We are dealing with relatively large family data (> 1,000 trios) with genotyping array data.
We would like to see the allele effects on children separately of mothers, fathers, and children, as in a previous study (PMID: 34282336).
Unfortunately, this is difficult for us to replicate, so we need another method.
Is there any software that can achieve this or a similar traditional approach?
Also, can that be applied to grandparents?
In my understanding of the above mentioned study is they:

identified the parental genomes within the child genome;
thereon conducted a GWAS for birth weight

Overall, the aim of the previous study was to estimate the alleles responsible for birth weight or in our situation a comparable phenotype. It is important to note the prior study did not use publicly available algorithms.
Thus, we wish to identify the maternal and paternal alleles within their progeny via public domain algorithms, and therefore do not necessarily need to exactly reproduce the previous study.
We already have imputed data, and are using GCTA for usual GWAS. To achieve our goal, we need to identify the transmitted and non-transmitted alleles in parents.

Hi @PenguinPartyH0, thanks for the extra info. Just to be clear I don't do eukaryotic genomics in general, but I've edited the question with your added information. It would be useful to explain what commercial algorithms the authors used, because it is likely someone here will know the freeware version. That is very likely the most important information required. You can include this info by editing the question. I personally don't think this is a problematic situation. I have never personally known any commercial algorithms that are more powerful than those in the public domain

Has your data been phased as well as imputed, or just imputed?

We have just imputed data and genotyped data.

The first haplotyping step is typically done in two stages:

Imputation of missing genotypes - can be done using the Michigan Imputation Server
Phasing of called genotypes to create semi-contiguous haplotypes - can be done using SHAPEIT

I have many thoughts about the best approach for GWAS (mostly around the idea of population subsampling), which could potentially also apply to what you want to do. Unfortunately, it's been a long time since I did family population genetics, and I've forgotten the things we did to deal with complex population structure. I'm pretty sure we used PLINK for at least some of it, but I notice that PLINK has a disclaimer about close family structure, which suggests it may not work so well for trio data:

This method does not properly adjust for small-scale family structure. As a consequence, it is usually necessary to prune close relations with e.g. --king-cutoff before using --glm for genome-wide association analysis. (Note that biobank data usually comes with a relationship-pruned sample ID list; you can use --keep on that list, instead of performing your own expensive --king-cutoff run.) If this throws out more samples than you'd like, consider using mixed model association software such as SAIGE, BOLT-LMM, GCTA, or FaST-LMM instead; or regenie's whole genome regression.


text:
BAF and LRR calculation and usage
I'm  working with a vcf file, in particular I'm trying to identify a gene duplication in a sequenced genome.
I studied bcftools cnv command, but for this tool is needed a vcf file containing a BAF and LRR values, that there are not present in my vcf file.
And so the first question is, can I calculate them starting from read depth and allele depth? Once that I have calculate BAF and LRR, can I modify the original vcf file inserting the calculated values? To have a vcf file containing BAF and LRR values I have to ask to the geneticist to use particular pipeline in NGS instrument?

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.


text:
PL and QUAL values on VCF file?
I want to filter my VCF file to include only the relevant information I need but I have some questions about the results I got.
This is what the first 10 lines of my VCF file looks like:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  DATA_VARIANT
chrI    826 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:344,30,0:57:10:0,10,0,0:0,5
chrI    831 .   C   T   71  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:139,36,274:71:12:0,8,0,4:3,2
chrI    836 .   G   T   59  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:140,49,373:59:16:0,0,12,4:4,1
chrI    848 .   C   T   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:556,48,0:75:16:0,0,0,16:0,5
chrI    855 .   G   A   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:399,48,139:59:16:12,0,4,0:1,4
chrI    858 .   ATT AT  255 .   CNV=1;TYPE=INDEL    GT:PL:GQ:DP:ADP:ACN 0/1:480,48,160:85:16:4,12:1,4
chrI    876 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:865,75,0:102:25:0,25,0,0:0,5
chrI    904 .   G   A   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:1261,112,0:138:37:37,0,0,0:0,5
chrI    920 .   A   T   67  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:244,145,1426:67:48:41,0,0,7:4,1
chrI    925 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:1786,157,1:255:52:0,52,0,0:0,5

QUAL:
As I understand it, the QUAL is based off of Phred quality scores, which would mean that a QUAL value of anything greater than 20 would be fairly sufficient in terms of quality for that genotype. However, most of my QUAL scores are 255 and from what I've read that means there is no mapping quality available?

What does this mean for my data? I used the same tool (NGSEP) to map and call the variants on this strain so it should have mapping quality scores.
Do I need to filter out all the 255 scores? Or can they still be used?

PL: From what I've read about PL scores, they indicate the likelihood of the presence of a genotype (0/0, 0/1, 1/1). Based on the results, there should be a clear least likely genotype, "almost" likely genotype, and a most likely genotype indicated by a value of 0.

Some of my PL values do not include a 0 (i.e., the second variant call = 139,36,274). What does this mean for my data?

Any help would be great! Thank you!
Note: I used NGSEP as my variant caller

The VCF specification provides the definition for the QUAL field.
However, QUAL values are often capped by variant callers to a given value. I have seen 100 being used and according to this documentation on NGSEP it seems 255 has been chosen for this caller as the maximum value. So you are all good keeping these high QUAL variant sites and filtering only those below a certain value (e.g. 20)
PL is also defined in the VCF specification: phred-scaled genotype likelihoods rounded to the closest integer (and otherwise defined precisely as
the GL field).
Again, different variant callers may have different implementations for the specification. GATK for example normalizes the values so that the most likely genotype has a PL of 0. But this might not be followed in other callers. It seems that NGSEP simply does not apply this normalization or does it differently. I think this is a detail you can ask directly the authors here by opening an "issue".
My suggestion is that you rather use GQ (genotype quality) if you want to hard-filter your variants. This is more standard to do, as opposed to filtering using PLs. These values are also in phred-scale.


text:
Filtering file with AWK and writing output to new file
Hi I'm trying to filter a column in a PAF file with the AWK command and write it to a new file:
awk '$12>=20 {print $0}' /PATH/FILE.paf > /OUT/PATH/FITLERED_FILE.paf
This command makes sure the 12th column has a value of at least 20 and then prints the entire contents of the file - with the filtered column - to a new file.
This works for one file but I have 90+ I need to filter through. I tried using a for loop:
for i in /PATH/TO/FILES/*_mapped.paf; do awk '$12>=20 {print $0 > "awk_{}.out"}' $i; done
Which didn't work. And then I tried using a pipe command:
cat /Users/rimo/Desktop/Pat_Assemblies/flye_mapped.txt | parallel "awk '{$12>=20 print $0}' /Users/rimo/Desktop/Pat_Assemblies/Assemblies_Flye/{} > awk_{}"
Which also didn't work. I think the problem I am having is format-based but I would appreciate it if someone could steer me in the right direction!

What operating system are you using? Do you have access to GNU awk (gawk on some systems, just awk on most Linuxes except recent Ubuntu and maybe some others).

I'm using a mac OS terminal

OK. Do you have gawk? No need to install it if you don't, it would just make some things easier if you already have it.

I don't have gawk but I could install it if it's easier! What would you suggest as an alternative?

You don't really need gawk here. I was mistakenly thinking you wanted to edit the original file which is made easier by gawk (you can use gawk -i inplace) and gawk won't choke if you have too many open files, but all you need to avoid that is to close the file handle after writing (see the last command in my answer), so you really don't need to worry. Just use the default awk.

The {} isn't in any way special to the shell or awk. It is used as a replacement by certain other programs, off the top of my head I can only think of find and xargs where it is used by convention, but it isn't a shell thing or an awk thing. Here, you are using it inside the awk command so it just breaks the awk's format because { } introduces a command block and this one is empty.
Your file name is $i, so if you want to make foo_mapped.paf into awk_foo_mapped.paf.out, you need awk_$i.out. Something like this:
for file in /PATH/TO/FILES/*_mapped.paf; do
    fileName="${file##*/}"
    awk '$12>=20' "$file" > awk_"$fileName".out
done

Note how I removed the print $0. This is because the default action for awk when something evaluates to true is to print the current line, so awk '$12>=20 {print $0 } is equivalent to awk '$12>=20.
You can also use your original approach of having awk deal with the file name by using awk's internal FILENAME variable. Like this:
awk '$12>=20{ print > "awk_"FILENAME".out" }' /PATH/TO/FILES/*_mapped.paf

Depending on the number of files, you might need to close the file handle after writing. If you get an error about too many open files, try:
mawk '$12>=20{print > "awk_"FILENAME".out"; close("awk_"FILENAME".out")}' /PATH/TO/FILES/*_mapped.paf


text:
Add file name to FASTA file header using SeqIO
I have several fasta files that represent the proteins from different orthogroups. However, the orthogroup that each of the proteins belongs to is not indicated within the files. I am leveraging BioPython and I want to add this information to the header of each sequence in the fasta files.
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
    for file in os.listdir(dirname):
        seq = SeqIO.parse(file, 'fasta')
        ortog = Path(file).stem
        for record in seq:
            record.description = ortog
        SeqIO.write(record, output, 'fasta')

However, the output always appears with the last sequence duplicated. In the case of three fasta files, the output appears as:
>BBC_8_672 OG0015396-PROT
MKFSSVTAITLATVATVATAKKGEHDFTTTLTLSSDGSLTTTTSTHTTHKYGKFNKTSKS
KTPNHTGTHKYGKFNKTSKSKTP
>CAT_6_3128 OG0020551-PROT
MEKIKEKLNSLKLESESWQEKYEELREQLKELEQSNTEKENEIKSLSAKNEQLDSEVEKL
ESQLSDTKQLAEDSNNLRSNNENYTKKNQDLEQQLEDSEAKLKEAMDKLKEADLNSEQMG
RRIVALEEERDEWEKKCEEFQSKYEEAQKELDEIANSLENL
>ZSA_0_3121 OG0025647-PROT
LLLQYQLHPALLPKALLLQYQLHQALLLKALLLQ*PALLLKALLLQYQLHQALLLKALLL
QYQLHPALLLKALLLQYQLHP
>ZSA_0_3121 OG0025647-PROT
LLLQYQLHPALLPKALLLQYQLHQALLLKALLLQ*PALLLKALLLQYQLHQALLLKALLL
QYQLHPALLLKALLLQYQLHP

Could anyone tell me why this is happening?

So what do you see when you add print(file) in the for file loop? Are some of the file names repeated?

This is because you are also parsing your output file. Remember, the first rule of debugging is to print all the things! Most bugs are because some variable doesn't actually have the value you think it has, and printing them is the best way of finding the problem. I reproduced what (I am guessing) is your setup here:
$ ls
OG0015396-PROT.fa  OG0020551-PROT.fa  OG0025647-PROT.fa

I then ran your script with an extra print statement:
#!/usr/bin/env python3
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
       for file in os.listdir("."):
              print("FF:",file)
              seq = SeqIO.parse(file, 'fasta')
              ortog = Path(file).stem
              for record in seq:
                     record.description = ortog
              SeqIO.write(record, output, 'fasta')

Which printed:
$ foo.py
FF: OG0015396-PROT.fa
FF: OG0020551-PROT.fa
FF: output.fasta
FF: OG0025647-PROT.fa

As you can see, since the output.fasta file is created by the open() call, by the time you use listdir, it's already there and included. This results in duplicated sequences (not always the last one). I'm sure there are more elegant solutions, but a simple workaround is to skip the output file name:
#!/usr/bin/env python3
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
       for file in os.listdir("."):
              if file == "output.fasta":
                     continue
              seq = SeqIO.parse(file, 'fasta')
              ortog = Path(file).stem
              for record in seq:
                     record.description = ortog
              SeqIO.write(record, output, 'fasta')

Of course, this whole thing is a bit overkill for such a simple task and you could just use sed and the shell instead:
for file in *fa; do sed "1s|$| ${file//.fa}|" "$file"; done > output.fasta

Thank you! It makes so much sense now. I managed to solve the problem by changing the identation of the SeqIO.write() function to the last for loop, but I didn't understand what was wrong before. Your resolution makes it clearer.

You're welcome @EduardoMenotti. But please also take M__'s advice, I'm sure he's right about using the id and he makes a very good point about not using seq as a variable name. Just make sure to skip the output file or simply save the output to a different directory or just write to standard output and redirect to a file.


text:
Understanding PAF format
Does anyone have any experience with PAF files (minimap2 output)? I don’t fully understand what my output means.
I’ve looked at the minimap2 explanation (https://github.com/lh3/miniasm/blob/master/PAF.md) but I want to use the data here to go back and extract regions in the genome that the gene target mapped to. I just don't know what specific columns to look at to accomplish this.
Here is a little snippet of a PAF file I have:
contig_25   522460  514241  515712  +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_26   413461  403718  406719  +   GENE_TARGET 4988    17      3649    1762  3633  60
contig_44   24218   16005   17476   +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_52   244411  236199  237670  +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_57   1043583 1025163 1026895 +   GENE_TARGET 4988    3046    4986    1376  2059  60

Is this saying that the gene target mapped to coordinates 514241 through 515712 on contig_25 but only coordinates 1962 through 3649 of the gene actually mapped to that contig? I guess I’m mostly confused on the descriptions for the coordinates and what they actually mean.
I think for my use case I have to use the gene target as the reference sequence as I am trying to extract homologous regions to sequences of interest. I.e., my project involves 90+ assemblies and I would like to extract the regions in the contigs where they map to the gene target. I used minimap2 to map the assemblies to the gene target, used sort and awk to get the best/high quality matches, converted PAF to BED format, and then used seqtk subseq to extract the sequences from the respective assemblies.
This is what I get when I use the gene target as a reference sequence:
GENE    4988    17      3649    +   contig_26   413461  403718  406719  1762    3633    60
...
GENE    4988    3046    4986    +   contig_57   1043583 1025163 1026895 1376    2059    60

Would I use columns 8 and 9 as the region range to go back and extract the sequences from?

As @user172818 mentioned in a comment to your question, it looks like you might have used the query as the reference; they should be switched round.
A minimap2 mapping command is typically represented as follows:
minimap2 reference.fa query.fa > results.paf

This follows the "indexed reference first" command line order used in the overwhelming majority of mapping tools. In the minimap2 PAF output (which you have linked), the first column represents the query sequence name, and the 6th column represents the target sequence name (i.e. the opposite order from the command line order, also following with established convention).
Based on the desired mapping operation that you have described (i.e. "I want to use the data here to go back and extract regions in the genome that the gene target mapped to"), your minimap2 command should probably look like this:
minimap2 genome_contigs.fa GENE_TARGET.fa > results.paf

This will appropriately map the target region to the regions in the genome, and then columns 8 and 9 will represent "regions in the genome that the gene target mapped to."
Minimap2 (and other mappings programs) match query sequences to a reference, with the usual expectation being that the output will contain a single primary mapping per input query sequence (sometimes with secondary mappings to other locations). In other words, they try to find one match to any reference sequence for each query. Changing the order of mapping to target-as-reference means that the program will attempt to find one match to the target sequence(s) for each input reference sequence, which is usually not what is desired.
In any case, the columns in the PAF file correspond to the closest sequence ID to the left of them (which is, confusingly, ordered the opposite way round to the command line order):
$ minimap2 [genome_contigs.fa] {GENE_TARGET.fa}
{GENE 4988 17   3649} + [contig_26 413461  403718  406719 ] 1762 3633 60
...
{GENE 4988 3046 4986} + [contig_57 1043583 1025163 1026895] 1376 2059 60

$ minimap2 [GENE_TARGET.fa] {genome_contigs.fa}
{contig_26 413461  403718  406719 } + [GENE 4988 17   3649] 1762 3633 60
...
{contig_57 1043583 1025163 1026895} + [GENE 4988 3046 4986] 1376 2059 60

that's very helpful! thank you!


text:
Is picking variants belonging to a specific group of genomes from VCF file possible?
GnomAD has publicly available VCF files that have variant data gathered from over 15k individuals. For a university project, I want to randomly select 1k of these individuals and get the variant data belgoning to them only.
I am very new to bioinformatics in general and upon checking VCF files and how they work, my current conclusion is that this task is not achievable through filtering the vcf file in any way since variants in these files have no direct connection to the genomes they were found in (no array of individual IDs in a variant entry that show which individuals had this variant).
But this task was given to me by a professor whom I believe is competent to realize the task was impossible if it is like I thought, so is there something I am missing?
How can I go about randomly selecting a 1000 individuals from this 15 000 pool, and then generating a new VCF file with the data from these selected 1000?
Are there any fields that show which individuals genome samples showed a spesific variant entry in the VCF file?

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.

When you say you want "1k of these genomes", I think you mean "you want the variants identified in 1k of these samples" and not the actual genomes. Is that correct? And what do you mean by "variants in these files have no direct connection to the genomes"? If this is a multi-sample VCF file, the zygosity field will tell you if this variant was found in a given sample.

Do you want to simply randomly sample 1K individuals from the VCF and produce a second VCF with the variant data of those 1K individuals ? That's what I understand,

@JRodrigoF Yes that is exactly what I want, it seems to me that I understood some spesifications of the vcf file wrong. Is this task doable, how would one go about achieving this?

It is. A VCF ought to be sub-settable to a specific list of samples. The only thing you have to come up with for your specific case, is a list of 1K random sample IDs within the original VCF.
Geenerally speaking, knowing that you can do this and how involves getting familiar with the VCF specification as well as the main bioinfo tool used to manipulate and work with this format: BCFtools
Let's use 1000 genomes data as an example.
You can use a bit of UNIX together with BCFtools in order to get a list of 1K random IDs from individuals:
(sort -R instead of shuf might work too)
bcftools query -l ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | shuf | head -1000 > myRandomIDs.txt

Then again BCFtools to subset the list of individuals from the original VCF and into a new VCF file:
bcftools view --samples-file myRandomIDs.txt ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz -o myNewVCF.vcf

If you dig a bit deeper into the BCFtools documentation, you will find many things you can do extra if you wanted, like re-calculating the allele frequencies based on this new cohort of individuals, if you wanted that too. Or subset for positions/variants, instead of individuals.
GnomAD data
Using GnomAD data, selecting 1K random samples is impossible since data is aggregated and instead of genotypes per individual (and individuals you can select thus) we only have the frequency of the variants across various sub-cohorts within the GnomAD project. No individual level data.
The most you can do is select one of such sub-cohorts and filter for variants with a frequency > 0 (or a range) in such cohort. You could get variants present only and only is a specific cohort too, by using either frequency information or other info such as allele counts, etc.
The example below filters for variants with a frequency above 0 in the South Asian cohort and 0 frequency in the Southern European cohort. You can have a look at the header of the VCF in order to get to know all cohorts included and the meaning of all the information contained within the INFO filed in the VCF
bcftools filter -i "AF_sas > 0 & AF_nfe_seu = 0" gnomad.exomes.r2.1.1.sites.21.vcf.bgz

But I'm afraid that's not exactly what you wanted to do in the first place.

The problemz is, while this works with the dataset you've given as the example, when I execute the bcftools query -l gnomad.genomes.r2.1.sites.chr21.vcf.gz it does not return any output. Does this indicate that in this VCF file does not include any sample ID information?

@ArdaCüce gnomad doesn’t actually contain individual level samples, it is a vcf of aggregated samples which gives you the allele frequencies (I.e. it contains averages), so it contains no sample level information. You might want to try a different dataset which does have individual level data, depending on the question you want to ask.

OK. I see, I for a moment thought your GnomAD VCF had sample-level data, but this makes sense actually. Please see if my updated answer provides a better answer,


text:
How to get full headers with tabs as separators with segmasker?
I try to pass a .fasta file with headers where tabs are used as separators as input file to NCBI BLAST+ segmasker:
>A0A059VAR9    IPR022879   V-type ATP synthase regulatory subunit B/beta
Unfortunately, the segmasker output's headers are shortened to the first tab:
>A0A059VAR9
Is there a way to get whole headers in the output file (without replacing all the tabs in the input file with white spaces)?

Could you supply a small data set please? You can fill is with dumbie amino acids ... I'll look at NCBI'ssegmasker local version because thats what I've got installed.

I think I also use local version. Here's a sample data: link

Thanks for the file, answered below. Sorry about the delay. This question isn't simplistic.

I used a local version of segmasker *. One way to solve this is,
perl -p -i -e 's/ //g' bioinfo.fa
perl -p -e 's/\t/^/g' bioinfo.fa | segmasker | perl -p -e 's/^/\t/g'  > bioinfo2.sm

Edit: I changed , for ^ because the OPs code already has , and will insert an extra tab (\t). Same change in the line below
Here the input file you presented in the link is called bioinfo.fa. If you need the pipe explaining let me know, but it's
fairly straight forward.
If you want to preserve space characters, i.e. ' ' or "\s" - because in the code above I delete them - then (not tested):
perl -p -e 's/\t/^/g' bioinfo.fa | perl -p -e 's/ /_/g' | segmasker | perl -p -e 's/_/ /g' | perl -p -e 's/^/\t/g'  > bioinfo2.sm

The issue the OP encounter is not unique to tabs (\t) but also occurs with spaces (\s), at least with makebasedb.
input
>A0A000 IPR004839   Aminotransferase, class I/classII   PF00155 41  381
GVSAIDGMNTLVRGQKLPIFSSSGQPHSELAMQIARQASVPEEEEGGDDEEGSEFAVIFGAMGITQEEANEFMEDFERTGALERSVVFMNLADDPAVERTVTPRMVLTTAEYLAFEKDYHVLVILTDMTNYCEALREIGAAREEVPGRRGYPG
>A0A000 IPR010961   Tetrapyrrole biosynthesis, 5-aminolevulinic acid synthase   TIGR01821   12  391
...

output
>A0A000 IPR004839       Aminotransferase        class_I/classII        PF00155 41      381
41 - 53
>A0A000 IPR010961       Tetrapyrrole_biosynthesis       _5aminolevulinic_acid_synthase  TIGR01821       12      391
41 - 53
>A0A000 IPR015421       Pyridoxal phosphate dependent transferase        major domainG3DSA:3.40.640.10  48      288
106 - 126
...

Please don't ask about solving this via an upstream makeblastdb because that gets complicated. The rate limiting step, i.e. why it is ignoring tabs is because of makeblastdb. Thats a consequence of the way its written ... its a developer issue (see blurb)

* Local segmasker installation ...

Install conda via Anaconda or miniconda
conda create -n blast3 python=3.11 (creating an environment ... yep 3.11 is out (wow!, if you don't understand why WOW! ask in chat)
conda activate blast3 (this is the conda environment, you can call it anything)
conda install -c bioconda blast

perl | segmasker | perl will now work. The Python installation is just good practice.

blurb about NCBI's blast The issue is makeblastdb, what happens is that long headers are a complete pain in algorithms because they use loads of RAM and you've got to carry them through reiterative loops (convergence etc ...). The way to solve it is via recoding (thats how I've done it). Essentially, you have a mini database forcing header recoding which activates at the start of the  algorithm and reassembles the original header at the end of the algorithm. The user never sees this step, but it's really important for efficiency. Thus within the reiteration the header isn't consuming RAM, because it's an internal code, beyond base 10 to squash the overhead even further. The thing is header can become a lot bigger than the sequence data and thats not cool when you've loads and loads of data, e.g. a blast database. NCBI just haven't bothered, like if you can't see it why do it attitude. What they've said is that 50 characters is the limitation of the header for makeblastdb... what NCBI is looking for is a unique id at the start of the header it then truncates everything beyond that, i.e. tabs, spaces ... anything it can grab. NCBI have carried this through into segmasker it's not needed, but makeblastdb and segmasker would normally be used together - thats why. In the OPs case this results in duplicate headers, so the file format falls over. This isn't the OPs fault.

Thank you for the exhaustive answer!

I agree, Python 3.11 has many great features. What's your biggest WOW there?

@maciejwww The criticism of Python is its slow. 3.11 is the major upgrade to focus on speed - its the fight back and may there be many more. Answer, you are welcome, note was a slight bug, I switched , for ^.


text:
How can I clean a PDB file using VMD or NAMD?
I have 3884 PDB files generated by pd2_ca2main. PD2 CA2main is server for reconstructing the backbone of Calpha proteins using Gaussian mixture models.
All of these files have duplicate atoms or alternate coordinates.
How can I clean these PDB files using Visual Molecular Dynamics (VMD) or Scaleable Molecular Dynamics (NAMD) *?
Can anyone tell me what command I should use in console mode?
* note, this server may be unavailable

You can use Chimera for this. The following code cleans your pdb file and write new pdb file.
import chimera
from DockPrep import prep
from chimera import runCommand, openModels

models = chimera.openModels.list(modelTypes=[chimera.Molecule])
prep(models)
runCommand("write format pdb 0 molecule.pdb")

You can run the code with
 chimera --nogui file.pdb code.py

how do I get chimera ?  pip install pychimera https://pypi.org/project/pychimera/ ??? is there a python3 version ?

https://www.cgl.ucsf.edu/chimera/download.html

https://pychimera.readthedocs.io/en/latest/ --> Use UCSF Chimera packages in any Python 2.7 interpreter


text:
Create VIP plot using R from the SIMCA VIP values
I am using SIMCA software to perform multivariate analysis and want a VIP scores plot in R as shown below.
which package is used for this type of plot (image retrieved from Google-images)??

Just a minutes VIP for SIMCO is Variable Influence Projection and the OP i assume wants to reproduce that in R because it is free. SIMCO is comercial.VIP has been used to quantufy transcriptomics via ML, but thats completely different VIP

Anyway the ML framework for VIP Variable Importance Plots in R is here https://cran.r-project.org/web/packages/vip/readme/README.html


text:
Good source for specific Isoforms overexpression in tumor cells (specifically, melanoma)
My group is interested in genetically engineering mice with label for the Slc45a2 gene, since "the protein is expressed in a high percentage of melanoma cell lines". (from Refseq)
How can I choose with high confidence the best isoform to use, i.e., the one that is most upregulated in melanoma cell lines?
Is there a curated database/ web-source of any kind that I can look it up?

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.

MelGene database, https://pubmed.ncbi.nlm.nih.gov/25407435/
and then more recent studies, e.g.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7927184/ . Just to note cancer isn't my thing, I just did a pubmed search.

Many Thanks. turns out this MelGene database no longer exists. at least I couldn't find any trace of it =\


text:
finding overlapping motifs to increase length of motif
I am able to find matching motif in my sequence, and I would like to now find overlapping motifs. Basically, after matching my motif, I want to find the 6 amino acids after it. This is the code below that I used to find the motif:
import Bio
import regex

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     result=regex.finditer(r"[YFWLIMVA]..[LMALVN]..[AGSTCD].[LAIVNFYMW]",sequence)
     for x in result:
     print(name, x.start(), x.end(), x.group())

The above code works perfectly becasue it give me the sequence id, positions and the motif. The output is below
P1  33 41 VTLLPAADL

Right now, what I want to do is to also get the 6 amino acids after matching this motif, such that I get an output like the one below where the length will now be 15,.
P1 33 47 VTLLPAADLLMAIID

The code that I have tried to get the 6 amino acids after my match is below.
import Bio
import regex

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     result=regex.finditer(r"[YFWLIMVA]..[LMALVN]..[AGSTCD].[LAIVNFYMW]",sequence)
     for x in result:
     print(name, x.start(), x.end() + 6, x.group())

This the output it gives me:
#It does not extend my motif by 6 amino acids, after getting the match, it only extends the
#end position by 6. From 41 to 47.
P1  33 47 VTLLPAADL

#My desired output is this which include the overlapping LMAIID, the LMAIID is after the
#match.
P1   33 47 VTLLPAADLLMAIID

I also tried the code below, but it returns an error.
import Bio
import regex

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     result=regex.finditer(r"[YFWLIMVA]..[LMALVN]..[AGSTCD].[LAIVNFYMW]",sequence)
     for x in result:
     print(name, x.start(), x.end() + 6, x.group() +6)

Cross-post https://www.biostars.org/p/9545486/#9545494

This is how it managed to give the desired output.
import Bio
import re

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     matches=re.finditer(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]",sequence)
     for x in matches:
     print(name, x.start(), x.end()+6,sequence[x.start():x.end()+6])

#Output
P1 33 47 VTLLPAADLLMAIID

This looks like a great answer, which has worked for you. Can you please flag this as your accepted answer (it's okay to do this for your own answer), so that the StackExchange automated algorithms do the right thing?

Just add 6 extra aminoacids to your motif:
matches=re.finditer(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]......",sequence)

That's all you need, now the existing code will give you the answer you want.


text:
Multiple values for memory resources shown when Snakemake is run
I am running a Snakemake (V7.14.0) pipeline on the cluster and provided resources in a resource config file like:
threads:
    tool1: 4
mem_gb:
    tool1: 12

Point to it in my rules like:
resources:
        mem_gb=config['mem_gb']['tool1']
threads:
        config['threads']['tool1']

And the snakemake cluster command like:
snakemake --configfiles \
    config/config_resources.yml \
    config/config_parameters.yml \
    --latency-wait 60 \
    --use-conda \
    --use-singularity \
    --printshellcmds \
    --jobname MyPipe{name}_{jobid} \
    --jobs 300 \
    --cluster \
    "bsub -W 200 -n {threads} -M {resources.mem_gb}G \
    -R 'rusage[mem={resources.mem_gb}G]' \
    -o run_logs/{name}_%J.out \
    -e run_logs/{name}_%J.err"

What I do not understand when I run the pipeline is that it shows this in the terminal:

threads: 4
resources: mem_mb=19442, disk_mb=19442, tmpdir=/tmp, mem_gb=12

So next to the one I provided (mem_gb), there are: mem_mb, and disk_mb. They are clearly not the same value. In the cluster job it looks as though the MEMLIMIT is indeed the one I give, but the MEM usage is different (higher). Is it because Snakemake does not understand the mem_gb and I specifically need to use that word?
Thanks in advance

You have defined mem_gb but not mem_mb in your resource configuration file. You can set them to same value by adding something like mem_mb=config['mem_gb']['tool1']. Additional information in their documentation.

Thanks for the response. So basically, Snakemake does not understand mem_gb. So to do this correctly, I need to actually *1000 the values. Because otherwise I get too low values (e.g. 6 mb instead of 6 gb).


text:
ETE4 manual or improved code description
ETE3 * (Environment for Tree Exploration) is a major package and in particular Python package to navigate phylogenetic trees. It is the premier coding library/package for tree manipulation.
ETE4.0.0 beta has been released via GitHub. The download contains an old ETE manual (2009) and a README file for ETE4, thats about it. Does anyone know the location of the new ETE4 manual, or list of functionality beyond what is in the README (posted below). Reading the code is the last resort ;-)

README ETE4

Read and write support for trees in Newick format
Multiple functions for traversing, searching and manipulating tree topologies and node annotations
Integration with NCBI Taxonomic database
Integration with GTDB database
Programmatic visualization framework
Visualization of huge trees based on adaptive zooming
Comparing trees
Phylogenetic functions:

orthology detection
phylogenetic distance?

Reference
* ETE 3: Reconstruction, analysis and visualization of phylogenomic data.
Jaime Huerta-Cepas, Francois Serra and Peer Bork.
Mol Biol Evol 2016; doi: 10.1093/molbev/msw046

The github page you link to says they have documentation on their official website (http://etetoolkit.org/) and, indeed, clicking on the documentation link there took me to  http://etetoolkit.org/docs/latest/index.html. Is that what you're looking for?

Thanks @terdon, this is the documentation for ETE Toolkit 3.0.0b34 (last updated Mar 18, 2016). I recognise the ETE3 to ETE4 is small outside phylogenetics, but it is a big upgrade. The authors may respond here, if not I'll send them this link.

I suspect you will get more traction if you just ask on their github. This isn't a question that is likely to be answerable by anyone other than the devs.

Good point. https://github.com/etetoolkit/ete/issues/651 ... will pop their response as an answer here. ETE4 has potential.

Just to feedback, the GitHub ticket was:

ETE4.0.0 beta contains an old ETE manual (2009) and a README file for ETE4. Is there an updated ETE4 manual beyond etetoolkit.org/docs/latest/index.html or the README of ETE4?

lm-jkominek wrote,

I know it's a bit of necro, but I wanted to second this request  - while the current docs work fine enough for ETE3, they have diverged quite far from ETE4 and there's a bunch of errors that come through when using the latest version that have no explanation...

Thus I'm not alone ... (nothing to do with UAPs BTW)
jhcepas wrote,

absolutely. We are finishing some final changes to the core of ETE4 and will start with the documentation. In the mean time, @jordibc has a list of most important changes regarding compatibility. Could you share somewhere as a wiki page @jordibc ?

jordibc replied,

Yes, we are working on an updated version of the manual to be released with ete4. In the meantime, most of the old code can be run with some of the changes mentioned in the part of the wiki on migrating from ete3 to ete4. @lm-jkominek and @sarscov2 thanks and hope this helps for the moment!

The GitHub wiki was updated here. The Wiki page is kept up to date,

"Jordi Burguet-Castell edited this page 4 days ago · 2 revisions"

The actual additional features (why ETE4 is needed) and how to access them are not listed here. However the changes in syntax are listed. Examples of changes ...
# ETE3
t = Tree(filename)
# ETE4
t = Tree(open(filename))
# ETE3
Tree(name='Hi', dist=1, support=1)
# ETE4
Tree({'name': 'Hi', 'dist': 1, 'support': 1})
# ETE3
Tree(nw, format=1)
# ETE4
Tree(nw, parser=1)

In summary there's different syntax. So we can start using ETE4


text:
How to resolve Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column while using merge function in R
I have two csv formatted data frames that I want to merge. First data frame contains 9 columns and the second contains three columns. Based on column names: start, end, function I am trying to merge the two data frames. Following is what I tried:
df1<-read.csv(file = "LuC_A_TE_gff.csv")
df2<-read.csv(file = "LuC_A_TE_1kb_range_ids.csv")
merged.data<-merge(df1, df2, by.x=c("start","end","function"), by.y=c("start","end","function"))

Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column.

df1

Output
chr,RM,sim,start,end,stat,strand,id,function
Chr01A,RepeatMasker,similarity,1,1122,4438,+,.,Target=rnd-4_family-4692 2 652;Class=LTR;PercDiv=4.9;PercDel=9.9;PercIns=5.4;
Chr01A,RepeatMasker,similarity,1123,2001,6668,+,.,Target=25:218931..219866_LTR 171 936;Class=LTR/Gypsy;PercDiv=6.5;PercDel=6.5;PercIns=0.7;
Chr01A,RepeatMasker,similarity,1998,2421,3164,+,.,Target=rnd-4_family-4692 102 559;Class=LTR;PercDiv=6.1;PercDel=8.0;PercIns=0.0;
Chr01A,RepeatMasker,similarity,2002,2119,1182,+,.,Target=rnd-5_family-5429 198 324;Class=DNA/Tc1-Mariner;PercDiv=5.1;PercDel=5.7;PercIns=0.0;
Chr01A,RepeatMasker,similarity,2284,3113,6265,+,.,Target=25:218931..219866_LTR 32 908;Class=LTR/Gypsy;PercDiv=4.6;PercDel=6.4;PercIns=0.7;
Chr01A,RepeatMasker,similarity,3003,3103,1238,+,.,Target=rnd-5_family-5429 159 202;Class=DNA/Tc1-Mariner;PercDiv=9.1;PercDel=5.5;PercIns=1.6;

.
df2

output
start,end,function
10000060,10000222,Target=rnd-6_family-1538 599 686;Class=DNA/CMC-EnSpm;PercDiv=24.8;PercDel=1.1;PercIns=3.2;
10000190,10000200,Target=33:34356884..34364107_INT-int 202 270;Class=LTR;PercDiv=20.5;PercDel=5.9;PercIns=1.6;
10000201,10000302,Target=rnd-6_family-3392 448 555;Class=DNA;PercDiv=10.9;PercDel=6.9;PercIns=0.9;
10000223,10000356,Target=rnd-6_family-1538 196 343;Class=DNA/CMC-EnSpm;PercDiv=20.3;PercDel=11.2;PercIns=0.7;
10000312,10000423,Target=fa.05_5_22233 5 115;Class=DNA;PercDiv=3.6;PercDel=0.0;PercIns=0.9;
10001256,10001332,Target=rnd-4_family-2291 200 273;Class=DNA;PercDiv=17.6;PercDel=0.0;PercIns=4.0;
10001590,10001909,Target=rnd-4_family-274 376 722;Class=DNA/CMC-EnSpm;PercDiv=12.3;PercDel=9.4;PercIns=0.9;
10001650,10001757,Target=rnd-6_family-7385 84 183;Class=DNA;PercDiv=16.3;PercDel=0.9;PercIns=9.0;
10002071,10002421,Target=rnd-3_family-11 116 430;Class=DNA/hAT-Ac;PercDiv=4.2;PercDel=0.0;PercIns=11.4;
10002933,10003171,Target=rnd-5_family-572 47 262;Class=DNA;PercDiv=17.4;PercDel=4.6;PercIns=15.7;
10003024,10003085,Target=Chr23:11562195..11564516 1491 1552;Class=RC/Helitron;PercDiv=4.8;PercDel=0.0;PercIns=0.0;
10003030,10003207,Target=rnd-4_family-825 65 227;Class=DNA;PercDiv=12.9;PercDel=7.3;PercIns=17.2;
10003132,10003142,Target=rnd-5_family-1284 903 950;Class=DNA/Kolobok;PercDiv=16.7;PercDel=7.7;PercIns=2.4;
10003143,10003941,Target=fa.24_3_1586 1 753;Class=DNA;PercDiv=17.3;PercDel=4.8;PercIns=11.2;
10003209,10003451,Target=rnd-5_family-502 491 715;Class=DNA/IS3EU;PercDiv=19.6;PercDel=7.4;PercIns=16.0;
10003271,10003496,Target=rnd-6_family-3401 557 818;Class=DNA/IS3EU;PercDiv=22.4;PercDel=19.5;PercIns=3.0;
10003414,10003567,Target=rnd-1_family-1158 5 166;Class=Unknown;PercDiv=25.1;PercDel=9.7;PercIns=4.3;
10003494,10003560,Target=rnd-4_family-825 2 63;Class=DNA;PercDiv=13.7;PercDel=4.5;PercIns=12.9;
10003631,10003854,Target=rnd-1_family-378 1 406;Class=DNA/P;PercDiv=4.2;PercDel=0.0;PercIns=0.1;
10004081,10004200,Target=fa.00_5_3944 125 245;Class=DNA;PercDiv=10.8;PercDel=0.8;PercIns=0.0;
10004147,10004282,Target=fa.13_3_23430 1409 1547;Class=DNA/hAT;PercDiv=6.7;PercDel=2.9;PercIns=0.7;
10004172,10004199,Target=fa.02_2_13330 434 524;Class=DNA;PercDiv=16.6;PercDel=3.1;PercIns=7.6;
10004200,10004427,Target=rnd-1_family-377 2 229;Class=DNA;PercDiv=11.3;PercDel=3.1;PercIns=3.1;
10004279,10004320,Target=fa.13_3_23430 882 923;Class=DNA/hAT;PercDiv=6.8;PercDel=0.0;PercIns=0.0;
10004321,10004656,Target=fa.15_1_6197 348 728;Class=DNA/Tc1-Mariner;PercDiv=13.1;PercDel=13.4;PercIns=0.0;
10004428,10004480,Target=fa.02_2_13330 525 709;Class=DNA;PercDiv=16.6;PercDel=3.1;PercIns=7.6;
10004482,10004804,Target=rnd-1_family-77 1 331;Class=SINE/tRNA-Core-RTE;PercDiv=7.5;PercDel=4.0;PercIns=1.5;
10004515,10004566,Target=rnd-6_family-2042 304 370;Class=DNA/TcMar-Tc1;PercDiv=7.7;PercDel=28.9;PercIns=0.0;
10004658,10005425,Target=rnd-5_family-3762 370 1163;Class=DNA/hAT-Ac;PercDiv=7.8;PercDel=3.8;PercIns=0.4;

When I loaded these files into R, I got the following result, which matches your error (except that it's been translated into German):
> df1 <- read.csv("df1.csv")
> df2 <- read.csv("df2.csv")
> merge(df1, df2, by=c("start","end","function"))
Fehler in fix.by(by.x, x) : 'by' muss gültige Spalte(n) spezifizieren

This error is obscure, and it took me a while to get my head around what was going on. The problem is that 'function' is a reserved word in R, so it renamed the function column to add a period [.] at the end of it:
> colnames(df1)
[1] "chr"       "RM"        "sim"       "start"     "end"       "stat"      "strand"
[8] "id"        "function."
> colnames(df2)
[1] "start"     "end"       "function."

So R was looking for a column name called "function", didn't find one, and got confused.
I think the easiest fix to this is to use the new names that R has given the tables:
> merge(df1, df2, by=c("start","end","function."))
[1] start     end       function. chr       RM        sim       stat      strand    id
<0 Zeilen> (oder row.names mit Länge 0)

[note: this has no output, because the examples you provided have no common start/end/function]
Alternatively, let R work out the common fields itself, as its default for merge is to intersect the column names (which is what you want in this case):
> merge(df1, df2)
[1] start     end       function. chr       RM        sim       stat      strand   id
<0 Zeilen> (oder row.names mit Länge 0)

Which, for completeness, is equivalent to this:
> merge(df1, df2, by=intersect(names(df1), names(df2)))
[1] start     end       function. chr       RM        sim       stat      strand    id
<0 Zeilen> (oder row.names mit Länge 0)


text:
How to identify mutations in a viral genome
I have a fasta file with multiple sequences comprising 38 sequences. The length of the sequences are around 11000 bp.
How can i get changes in the genomes based in a reference genome? (aa subtitutions or deletions)

Is this a 10700 genome by any chance?

Python3 via the difflib library/package. AlignIO can be leveraged (Biopython) to construct the align list. To be honest if your coding isn't strong you could just cut and paste the alignment into this code. Nucleotides not amino acids ... same idea just easier to type. It will work for aa and gaps. "Ref" (output) refers to the first alignment position only.
import difflib

align=['AGCTAGATGATATGA', 'AGCTAGATGATATGA', 'AGCTAGGTGATATGA', 'AGCTAGATGATATGA', 'AGCTAGATGATATGA', 'AGCTAGAAGATATGA', 'AGCTAGAGGATATGA', 'AGCTAGATCATATGT']
ref= 'AGCTAGATGATATGA'
for query in align:
    print('{} (ref) => {} (query)'.format(ref,query))
    for loc,pos in enumerate(difflib.ndiff(ref, query)):
        if pos[0]==' ':
            continue
        elif pos[0]=='-':
            print(u'Ref "{}"'.format(pos[-1]))
        elif pos[0]=='+':
            print(u'Query has "{}" at position {}'.format(pos[-1],loc))
    print ('')

Output,
AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGGTGATATGA (query)
Ref "A"
Query has "G" at position 7

AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGAAGATATGA (query)
Ref "T"
Query has "A" at position 8

AGCTAGATGATATGA (ref) => AGCTAGAGGATATGA (query)
Ref "T"
Query has "G" at position 8

AGCTAGATGATATGA (ref) => AGCTAGATCATATGA (query)
Ref "G"
Query has "C" at position 9
Query has "T" at position 16

Sequences must be aligned, i.e. ref sequence within the global alignment and separated. Hope that makes sense. It will work with multiple mutations between ref and query as the last example points out (two mutations per example query was just too much to type).
A reproducible example in future please.


text:
Lifting snp rsid build 37 to 38
I’m looking to lift over a dataset of GWAS summary stats from build 37 to 38, including converting rsid. Downstream it is to be used for co-localisation with another dataset of build 38.
Would anyone be able to help with code for doing this please?

Have you tried some tool or approach already?

Should we assume you are working on human data? Are the SNPs in coding regions only or do you also have intronic or, even worse, intergenic? Please [edit] your question and give us some context. Ideally also an example of your input.

There are multiple tools to use for this. Might I suggest the following reads: Introduction to liftover tools, where they also refer to this usefull post BioStars post. Maybe also consider reading this paper for a word of caution: Paper.
Maybe you can also clarify you question more, if these links do not help you with what you need.

I have recently worked on a couple of tools to do exactly this. They are +munge to convert summary statistics into files following the GWAS-VCF specification and +liftover to convert a VCF from one genome build to another. This will automatically take care of REF⇆ALT switches by flipping the effect size and it will correctly handle indels, if there are any, in your summary statistics file


text:
What options should I explore to improve the output of the protein reconstruction?
We developed a neural network-based protein reconstruction tool to reconstruct the main chain from only CA atoms.

we generated data from some selected PDBs from the RCSB website to train an NN model.
we then use those data to train the NN model.
we select some test PDBs, strip all atoms except CA atoms, and save them in files.
we pass those CA-only PDBs through the NN model, obtain a reconstructed main chain, and save them in files.
we compare original and reconstructed PDB files and calculate CRMSD values.

If we obtain a CRMSD value of 0.3559 and want to improve it what options should we explore?

Using Keras.
Not CNN nor RNN. It's a straightforward MLP (4 layers in total). No feedback, no convolution.
training set is actually large: 1398438 rows and 102 columns in the data file.

The first step is to assess the accuracy, recall and precision of the model then look at confusion matrix. Thereafter consider ROC AUC and if all those are fine ... it's biological. If the model is not 'accurate' where accuracy >0.85 (minimum) its about training the model and that is a different question.
My guess (edit: which from the comments was wrong) is:

the sample size is insufficient for a presumed deep learning model.
this will generate problems in the accuracy of the model.

The solution to low sample is called 'data augmentation', a lot of biological knowledge is needed to successfully augment a data set. In this case I suspect the level of biological understanding is strong, but the understanding of deep learning is less advanced.
In summary, low accuracy is a feature of small sample sizes within a deep learning model. To give some idea deep learning needs 500 000 sample size to be reasonable.

Update from the comments, the data set is large:

a sample of 1398438 rows and 102 columns (data points).

This is a good sample for deep learning and low sample size is not the cause.
The model is keras via a 4-layered MLP, so thats fine, but not clear whether its CNN (probably) or RNN.
The likely issue is under-fitting and this is diagnosed via the above statistics for example accuracy and ROC AUC and solution is to enhance the MLP. If there is over-fitting occurring thats a different question.


text:
Implement a WDL workflow for a docker container I built
I start saying my background is Biology and I'm totally new to WDL, so I'm aware my code will be messy but I'm trying to learn.
My problem is, essentially, to run a Docker container I built and tested for a tool on WDL (Terra platform). Now, I'm struggling a bit in understanding how to get input and generate output, and if they have to be contained in a path within the Docker. I successfully run this tool on a normal cluster, but Terra is a different story...
Shortly, what the tool does is determining the k-mers distribution in a pangenome and use the haplotypes support from the graph structure to infer genotypes for samples "mapped" on the pangenome graph space (more info here).
The tool requires as input a VCF file from the pangenome, a reference genome and a set of samples to genotype; all this information is specified within a .yaml file the user would normally edit once the tool is installed on a cluster to point to the correct directories where files are.

Firstly, I would like to give users the choice over what pangenome VCF file is loaded, as there are two main ways to build pangenome graphs and the genome inference is substantially affected by them.
Secondly, I want also to give a choice over the type of reference used, as variants for the samples can be called both against GRCh38 as well as CHM13.
Finally, I would like to give users freedom over the number of samples to genotype but also understand where this data (uncompressed .fastq files) have to be access on Terra and loaded in WDL

The tools generates five directories after the execution is completed, I'm not sure whether I have to create them beforehand for output to be access by users. Let me know.
Following, the script I wrote up to now. I have a few problems at different stages of the implementation as you can see. Please let me know if you can help and what kind of additional information you might need. Thanks in advance!
code version 1.0

workflow PanGenie {
    input {
        File? PANGENOME_VCF # input vcf with variants to be genotyped
        File? REFERENCE_GENOME # reference genome used to call variants from the graph
        Array[File]? INUPT_READS_SET # actual .fastq file
        Array[String]? SAMPLEs_FASTQ_READS # reads (FASTA/FASTQ format, uncompressed) for samples' to genotype

        Int CORES = 24 # number of cores to allocate for PanGenie execution
        Int DISK_STORAGE = 100 # storage memory for output files
    }

    File PANGENOME_VCF ## lets the user decide whether to use a mc-CACTUS or a PGGB based vcf file
    File REFERENCE_GENOME ## lets the user chose what reference to call variants from (GRCh38 or CHM13)
    File PANGENIE = "/app/pangenie/build/src/PanGenie" ## path to executable within Docker

    task pangenie_exe {
        input {
            File PANGENOME_VCF
            File REFERENCE_GENOME
            Array[File] INUPT_READS_SET
        }

        command {
            snakemake --cores CORES
        }

        output {
            Directory benchmark = ""
            Directory genotypes = ""
            Directory input_vcf = ""
            Directory pangenie = ""
            Directory pangenome = ""
        }

        runtime {
            memory: "300 GB"
            CPU: CORES
            disks: "local-disk " + DISK_STORAGE + " SSD"
            docker: "quay.io/overcraft90/eblerjana_pangenie"
            preemptible: 1 # can be useful for tools which execute sequential steps in a pipeline generating intermediate outputs
        }
    }
}


text:
conversion of a gbff file to a gff file
I want to get the annotations in gff format for a assembled genome. But for some species there is no available annotation in gff format in NCBI but only in gbff format. And I don’t know how to convert a gbff file to a gff file. So, I am seeking your help. Can I have your suggestions here? Your guidence will be appreciated.

answered here: https://www.biostars.org/p/9545716/#9545750

Thanks for your response.  Is the following link be applicable?
https://github.com/Tesla-huang/Scripts/blob/master/gbff2gff.pl

https://github.com/jorvis/biocode/blob/master/gff/convert_genbank_to_gff3.py

I don't know either way (just happened to spot GenoMax's answer at Biostars)... but if either/both of those links do the trick for you, then maybe that answer was wrong and this is possible.


text:
finding the frequency of a motif-matching peptide
This question has also been asked on Biostars
I want to find the frequency of each motif-matching peptide.
The sequences are:
> Sequence_1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence_2
MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR

>Sequence_3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

> Sequence_4
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence_5
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

The fasta file actually has more sequences, I just chose the top 5, to put on here.
The code used is below:
import Bio
import regex

from Bio import SeqIO
input_file = 'Sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     matches=regex.finditer(r"(P[A-Z]{2}P..)",sequence)
     for m in matches:
        print(name, m.start(), m.end(), m.group())

This is giving me my desired output, which is below:
Sequence1 73 79 PFLPGK
Sequence1 281 287 PLPPKL
Sequence1 288 294 PSSPAH
Sequence4 73 79 PFLPGK
Sequence4 281 287 PLPPKL
Sequence4 288 294 PSSPAH

Now I want to:

Calculate the number of sequences in that file. I know in linux environment I would have to use grep -c ">" Sequences.fasta to get the total number of sequences in the fasta file, but how do I do it in python?.

Once I know how many sequences are in file in total, I want to count how many of those sequences have the motif-matching peptides above.
For example how many sequences have PFLPGK at positions 73 79? In the example above PFLPGK its in sequence 1 & 4, so the frequency is 2, meaning that two sequences have PFLPGK. How do I find that frequency for each of the peptides using python.

Many thanks for your assistance.

Hi @thole the historic protocol on the site was to ask an OP what code has already been written to address in this instance Q 1 and 2. I personally think there should be a grace period. However, please avoid asking the site to write your code base from scratch for future reference.

Yes, pelase don't treat this site as a free script writing service. We are happy to help you with your work, but we won't just write it for you like this.

Just to add, this just for Q2. Q1 is a fair question by itself any time and I like the solution.

@M__ thank you, I was able to finally figure out how to work on Q1 without using Biopython. As for Q2 I am still working on it. Do you by any chance have any suggestions on books or free sites that I can use to continue learn python? Once again thank you for your suggestions

@thole post updated as requested, the information sort is below the ------. There was a bug in my implementation of Q1 it now works.

For the first part is easy,
from Bio import SeqIO
seqs = SeqIO.to_dict(SeqIO.parse("myfile.fa", "fasta"))
print (len(seqs))

The second part is solvable Python's amazing slicing, this is the concept ....
from Bio import SeqIO
loci = {'PFLPGK':[1,10], 'PFLPGK':[50:60] }
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
mydict = {}
for key, values in loci.items():
    myrange = '-'.join(values)
    mydict = {key:myrange} # you need to initiate the dictionary without defaultdict
    for record in fasta_sequences:
        if record.seq[values[0]-1:values[1]-1] == key:
             mydict[key][myrange].append(record.id)

To print it:
key, value in mydict.items():
    print("No. of" + key + '=' + len(mydict[key][value]))

Note 1, the second answer is concept only.
Note 2 there can be an issue between the list location and the sequence location ... that is super important otherwise it will not work. I've added -1 from both locations of the amino acid and the location on the list, but this needs checking. This is because lists start at 0 whilst amino acids start at 1

Update Advice on Python requested by OP

Answer to Q1 has been debugged, tested and works.
Use Q2 as the basis ... just debug it (see point 4.), it will be a bit fiddly because moving from lists to aa strings.
Use an IDE, I use Spyder: don't get too caught up in "the best IDE" they just make life easier.
Book: Mark Lutz "Learning Python 5th Edition" by O'Reilly section on lists (slicing) and dictionaries *.
Online resources GeekforGeeks on lists (slicing) and dictionaries.
BioPython docs on SeqIO search for record.id and record.seq ... you need to understand the SeqIO object.

* This is physically a huge book so an electronic version is strongly preferred just for logistics. This goes beyond the full modern implementation of OOP in bioinformatics (which is not very modern in truth and can become unreadable).


text:
LD block labels annotation for GWAS summary data
I have a txt file summarising the result of a GWAS on an European population. Its structure is the next one:
data = data.frame(chr = c('1', '1', '1', '1', '1'),
              bp = c(740098, 787889, 952073, 993492),
              snp = c('rs12138618',  'rs4951864',  'rs3128126',  'rs4075116'),
              p = c(0.1, 0.04, 0.7, 0.9))

Imagine the first two snp listed (rs12138618, rs4951864) belong to the same LD-block (for r^2 = 0.1 or similar criteria). Where can I get this information from (table/database) or what package in R can I use to label these dependent SNPs and obtain something as:
   chr    bp        snp    p    ld_block
1   1 740098 rs12138618 0.10        21
2   1 787889  rs4951864 0.04        21
3   1 952073  rs3128126 0.70        29
4   1 993492  rs4075116 0.90        34

I am using big data, so something computationally-wise is preferred.
Note: Obviously I am making up the numbers, but I hope my point is clear. I am also new to GWAS, so it is very likely that I can be missing/misunderstanding concepts.
What if I do not have the genotype data but I know that my population is CEU? Is there any reference-genotype data that I can provide to plink instead?

Is your overall goal to find the 'best' SNP within each block of LD, or do you explicitly need to know all SNPs within some LD block for another reason? Can you clarify what the overall goal is please.

I want to explicitly know to which LD block belongs each SNP, in order to obtain something similar to the represented output stated in my question. My problem is that I do not know from where can I get some "reference-genotype data" for an European population (I thought 1000G may have something but I did not manage to find it)

@Gero what if one SNP is in LD with more than one LD block? Most frequently SNPs are in continuous LD with other SNPs upstream and downstream, and you can't neatly split them, though there are of course inference techniques for segmentation. You can infer "haplotype blocks", which are one way of segmenting this information. Usually this is somewhat data dependent, do you want to do this according to your specific data using a tool or using some global reference (e.g. tag SNPs from HapMap)?

My plan was to: 1) Do a LD matrix, 2) find all SNP pairs based on r^2,  3) Create LD blocks considering LD as a transitive property (if one SNP is in LD with more than one block then I will consider those "2 blocks" as 1 big unique block). Even if LD is continuous as long as the distribution of the SNPs in the genome is not homogenous I guess many blocks will be obtained, hopefully ∼1,000,000

On the other hand, thank you for the tools. I'll try to investigate how they work and do the maths for the segmentation! (P.S. My plan is to use a global/population reference).


text:
Why do Illumina 850k/EPIC arrays ignore CpGs which are "GC" in the forward strand?
CpGs are symmetrical, in that a CG sequence on the forward strand is hybridized to a GC --- and both dinucleotides on each opposing strand are CpGs dinucleotides which can be methylated. Conversely, CpGs can be GC on the forward strand but CG on the reverse strand.
FORWARD -> 5'--CG--3'  [OR]  5'--GC--3' <- FORWARD
REVERSE -> 3'--GC--5'        3'--CG--5' <- REVERSE

The assignment of "forward" and "reverse" strandedness is more or less arbitrary.
Given the above, why does it seem like the Illumina 850k (aka EPIC) array only profiles methylation from CpGs which are CG in the forward strand, while ignoring CpGs which are GC in the forward strand? I would also love to hear if my premises are wrong.
Here is some R code illustrating the issue at hand:
suppressPackageStartupMessages({
  library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
  library(tidyverse)})
data(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

IlluminaHumanMethylationEPICanno.ilm10b4.hg19 %>%
  getAnnotation() %>%
  as_tibble() %>%
  count(forward_seq=str_extract(Forward_Sequence, "\\[[ATCG]{2}\\]"))

# Results:
# # A tibble: 3 × 2
#   forward_seq      n
#   <chr>        <int>
# 1 [CA]          2922
# 2 [CG]        862927
# 3 [CT]            10
```

I cross-posted this to Biostars: https://www.biostars.org/p/9545804/

My question was kindly answered by Biostars user Papyrus, who pointed out that a GC in the 5' to 3' direction is a GpC dinucleotide, not a CpG dinucleotide. GpC dinucleotides are then biologically and chemically distinct from CpG dinucleotides, due to the directionality of the DNA polymer.


text:
How to wget a gff from NCBI database?
Say I want to use wget to download a gff file for the X chromosome of release 5.57 of the Drosophila melanogaster genome. The NCBI page with the gff file for just this is here.
However, I do not see any way to directly wget the gff from this page. What I can do is click "Send to" at the top, select File, and change the Format dropdown to GFF3, and finally click 'Create File', and the gff file downloads onto my computer.
But then I lose the ability to use the -a flag to log the download details, I can't download the gff using wget from a remote server directly but instead I need to download it locally and then upload it, etc. Is there any way to use wget to download gff files from NCBI programmatically instead of having to download it manually from the website?

Are you looking for these?

@user172818 If you mean my example, not exactly, I'm looking for release a 5.57 gff based on chromosome X, that link lets me download a release 6 whole genome gff. But I think I figured out what to do, and I posted an answer to my question.

my suggestion would be to go via the assembly page.
For your example link, I would

Click the link to the BioProject: PRJNA164
On the right side, under "Related information" choose the assembly
On the assembly page, under "Access the data" on the right side, I would choose "FTP directory for RefSeq assembly"

Once in the RefSeq folder, copy the link to the *_genomic.gff.gz and use wget.
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/215/GCF_000001215.4_Release_6_plus_ISO1_MT/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.gff.gz

Edit: The initial question requested an older versions, which can be found in RefSeq's archive of prior releases. For this specific Drosophila assembly chromosome X, the version 5 can be found at https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5/CHR_X/NC_004354.gff

This does show how to wget from NCBI, although it wouldn't work on the specific case I mentioned. This gives me the whole genome, although I'm curious on getting a single chromosome. Also it will give me release 6 (the latest one), whereas my scenario was for release 5. However I think I found a way to do this.

Pardon, the RefSeq archive seems what you are looking for. Have a look at https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5/.  The path to the X chromosome's gff should be straight forward. https://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5/CHR_X/NC_004354.gff

That's exactly it! Could you edit these links into your main answer so I can mark your comment as accepted?

You're welcome. Glad to have helped

I believe I figured out a way to do this for my particular example (although in general cases Adrian's answer works). To get release 5.57 (an earlier release) of chromosome X of Drosophila melanogaster, I believe you can go to the Datasets Genome page of D. melanogaster, scroll down to the Revision History and select Release 5, and at the top of the page, NCBI lets you click on a 'datasets' or 'curl' box to copy the exact command line argument needed to download a GFF (and other files) of that whole genome. Then, I believe it should be possible to open the gff in Python using a package like gffutils, and filter only for chromosome X rows, and then write those rows to a new gff file.

If you are willing to filter, there is no need for python or writing a whole script. A simple grep -P '^NC_004354.3\t' foo.gff > chrX.gff will do.

Did not realize I could do that, thanks. Although you'd probably have to add a bit to this command to also preserve the header lines. Maybe (grep "^#" foo.gff; grep -P '^NC_004354.3\t' foo.gff) > chrX.gff


text:
Understanding regularization step in Meta-cell pre processing
I'm reading through this paper: Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., ... & Tanay, A. (2019). MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions. Genome biology, 20(1), 1-19.

Maybe I am not familiar with this notation but I don't understand the regularization procedure as follows:

What is happening in these max functions? What does the 0 indicate? They have some terms $x$ and then in the max function they write $max( x, 0)$.
I assume $K$ is some regularization parameter, like a cost. Is $s_{ij} * s_{ji}$ the rank product?  I.e. the geometric mean?

Without looking in too much detail:

We may suppose that $S$ is a matrix with dimensions indexed as $(i, j)$. So $s_{ij}$ is the cell of $S$ indexed to the $i$th row and $j$th column, and $s_{ji}$ is the cell vice versa. Therefore, symmetrizing $S$ is a matter of ensuring that $s_{ij} == s_{ji}$, which you can achieve by multiplying the (ranked?) entries. Symmetric matrices are nice for doing certain kinds of math, and yield a handy representation of undirected graphs, which the authors appear to be using for kNN.
$max(x, 0)$ is just saying that you want either $x$ if $x$ is non-negative, or $0$ if it is negative. Thus, you make the whole matrix non-negative, which makes some handy matrix math possible. You can see how this works computationally by inspecting e.g. the Python max() function.