text:
Generating 3D coordinates error
I have been trying to convert a multi entry sdf file to pdb using openbabel. Everytime I check the Generate 3D coordinates box I end up with .ENT files instead of .PDB files. What am I doing wrong.
Thank you.

An ent file is a pdb file. Just change the extension. PDBe for example provides their PDB files as such.

ok. I tried opening the .ent file with Autodock but I kept getting the error, .ent file is an unrecognized format. How can I simultaneously change the extension for 50,000 ent files to pdb without having to do it manually for each one individually.
Thank you,
Dr. Hussain.

In your file manager, select all and rename


text:
Creating custom VCF file
I'm stuck on what seems to be a banal problem, but I can't work out where the issue is. I need to create some customised VCF files that specify particular variants of interest.
This is a custom VCF I created for this variant:
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20240308
##source=Ensembl
##reference=GRCh38
#CHROM  POS ID  REF ALT QUAL  FILTER    INFO
20  4699752 rs74315403  G   A   .   .  .

The formatting seems correct to me; the values are separated by tabs, all relevant column headers are present, there is a value for every header. I try to use it in the Mutect2 --alleles and -L options, as follows:
gatk Mutect2 \
  -R ref.fa \
  --alleles my.vcf.gz  \
  -L my.vcf.gz\
  -I input.bam \
  -O output.vcf

This gives me the following error:
Caused by: htsjdk.tribble.TribbleException$MalformedFeatureFile: Unable to parse header with error: Your input file has a malformed header: there are not enough columns present in the header line: #CHROM     POS     ID      REF     ALT     QUAL  FILTER    INFO, for input source: path/to/my.vcf
        at htsjdk.tribble.TribbleIndexedFeatureReader.readHeader(TribbleIndexedFeatureReader.java:265)
        at htsjdk.tribble.TribbleIndexedFeatureReader.<init>(TribbleIndexedFeatureReader.java:104)
        at htsjdk.tribble.TribbleIndexedFeatureReader.<init>(TribbleIndexedFeatureReader.java:129)
        at htsjdk.tribble.AbstractFeatureReader.getFeatureReader(AbstractFeatureReader.java:121)
        at org.broadinstitute.hellbender.engine.FeatureDataSource.getTribbleFeatureReader(FeatureDataSource.java:433)
        ... 15 more
Caused by: htsjdk.tribble.TribbleException$InvalidHeader: Your input file has a malformed header: there are not enough columns present in the header line: #CHROM       POS     ID      REF    ALT      QUAL  FILTER    INFO
        at htsjdk.variant.vcf.AbstractVCFCodec.parseHeaderFromLines(AbstractVCFCodec.java:149)
        at htsjdk.variant.vcf.VCFCodec.readActualHeader(VCFCodec.java:117)
        at htsjdk.tribble.AsciiFeatureCodec.readHeader(AsciiFeatureCodec.java:79)
        at htsjdk.tribble.AsciiFeatureCodec.readHeader(AsciiFeatureCodec.java:37)
        at htsjdk.tribble.TribbleIndexedFeatureReader.readHeader(TribbleIndexedFeatureReader.java:263)

Supposedly there's something wrong with the header. Can you spot the issue? Thanks!

That sounds like your header isn't tab separated. You can check for the presence of tabs with od -c. So, to check the header, do
zgrep '^#C' my.vcf.gz | od -c

Then check the output, it should show \t between each header field, like this:
$ zgrep '^#C'  variants.vcf.gz | od -c
0000000   #   C   H   R   O   M  \t   P   O   S  \t   I   D  \t   R   E
0000020   F  \t   A   L   T  \t   Q   U   A   L  \t   F   I   L   T   E
0000040   R  \t   I   N   F   O  \t   F   O   R   M   A   T  \t   E   D
0000060   2   0   0   1   3   9   6  \n
0000070

Nice on the od -c. I generally use cat -te, but this looks easier to eyeball

source: Stack Exchange

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

sort: votes


text:
How do you find differences when your metagenome sample size is small?
I received a matrix of genomics features (KEGG annotated metagenomes) abundance from 6 samples belonging to 2 groups. My collaborator is interested in finding differences between the two groups.
The matrix, after scaling and removing zero-variance, is roughly 4000 features. Via shapiro test I already know some 30% of the features are not normally distributed across the samples, so I assume PCA is out of option and any test with the normal distribution assumption.
I already made some preliminary analysis using NMDS to cluster my samples and was thinking of trying also t-SNE and recursive feature elimination. However I am not sure how to find "significant" differences given the small sample size. Any suggestion? Would something like LEfSe work in this instance?

I think there are several ways to figure this out. I wouldn't worry too much about normality as most of the methods are either nonparametric or pretty robust to violations.

Since you are using KEGG, I would suggest working at the pathway or module level of annotation (if your features are e.g. KEGG Orthology groups). In principle, you can then combine information across features as well as across samples. Simplest I can think of, without even looking at the literature, would be to use permutations (shuffle row and column labels) to look for modules different in abundance between groups. (Many existing methods are doing some more sophisticated version of this.)
I think that LEfSe in theory should work quite well here, it's been a while since I thought about it but I believe this is the problem that it is designed to solve, assuming that I understand your question.
There are a lot of methods that are previously described for doing this kind of work. I would suggest reading some reviews of methods, for example the section "Multivariate differential abundance testing" in this article, or here, looking for some inspiration. To some extent it depends on which questions you are interested in answering and what exactly your features represent.

Calle's article has been greatly helpful! I did not know of all the work done on compositional data analysis. Now I have plenty of new tests to try. Thanks.

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


text:
Running GSEA or comparable analysis with multiple variable interaction?
I am working with an RNA-seq dataset generated from a complex experimental design. Our main question is focused on the interaction between two variables, one (INJ) with two conditions (L,S) and the other (SOC) with four conditions (A,B,C,D). We also have a third variable for tissue type (TIS) from 3 tissues (H,N,T), although we are not interested in this variable for this particular question.
We were able extract DE genes underlying the interaction of interest in DESeq2 using an LRT test. We also ran WGCNA to get gene clusters and looked at trait correlations. Using these WGNCA clusters and VST normalized count data, I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables.
What I've tried so far was generating a pseudo-variable that is a combination of INJ and SOC with 8 separate conditions for running in GSEA, but this doesn't quite get at the interaction between variables since we really want to see if the effect of INJ differs between SOC conditions.
One thing I was thinking was running a GSEA-like analysis on log2 ratios between INJ states for each SOC condition, but these would have to be summary statistics for groups of samples and I would lose the variance of expression within groups (maybe that doesn't matter though for exploratory analysis?). This would also reduce my "sample" count to 1 per condition, which isn't ideal for this type of analysis and I would probably have to look at differences between gene sets instead of phenotypes for more robust stats.
Any thoughts on ways to do this? I don't need to use GSEA if there is another analysis I could run on gene sets that allows for more complex experimental models. Thanks in advance.

You could read about differential co-expression. But that would be feasible depending on the number of samples you have for each condition...

I'm not sure I understand "I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables". I would run an association analysis (regression) for module eigengenes and the appropriate interaction term. The (most) significant modules are your candidates. This step is simply regression (in the simplest case a linear model should be sufficient), no GSEA needed. You can then subject the identified modules to GSEA or similar to get an idea of their biological function.

Thanks for the input and sorry for the late response. I've been looking into association analyses for modules/eigengenes and I'm thinking of running LRT analyses using GLMs with and without the interaction term of interest. Would you suggest I run separate analyses using each module's eigengene values individually (vs phenotypes) then correcting for multiple hypothesis, or can you recommend a more elegant way to test all modules at once? We would also like to run 3 separate analyses for each of 20 modules based on tissue sample subsets, so running 60 different analyses is not ideal. Thanks.

For module eigengenes, a linear model such as ones implemented in limma should be appropriate, no reason to run GLMs. With a linear model, you can simply test the significance of the interaction coefficient which is output by all linear model functions that I know (as long as the full and reduced model only differ by one coefficient). As far as I am aware, the only way to run such association analysis is fitting the models to the eigengenes one by one (perhaps with some empirical Bayes modifications as in limma) and then correcting for multiple comparisons.


text:
Decontaminating RNA seqs for de novo transcriptome assembly and annotation of novel eukaryotes
I have raw paired-end RNA-seq reads for two novel eukaryotic species. Some background: the reads represent a copepod (arthropod) species each. The mRNA for each read set was obtained by extracting mRNA from multiple individuals. I am not exactly sure what sequencer was used (some Illumina short read machine as far as I can tell; I am not directly involved in this project).
I'd like to assemble these reads into transcriptomes de novo as no reference genomes are available as far as I know; my hope is to assemble with Trinity and then annotate with Dammit using its default databases (and check for completeness w/ BUSCO).
However, it looks like something went wrong during the sequencing runs (or during sample preparation) because the sequence GC content plots from FastQC look like this:

Here the curve with the tall peak is one sample, and the more flatter curve is the other. There are two closely related curves per sample because FastQC apparently seems to have recognized each read (forward and reverse) as a separate sample.
I first presumed that this is just rRNA contamination (I believe the mRNA was prepared via poly-A selection), and I filtered the reads against all of sortmerna's default databases (which I believe include sequences from SILVA and Rfam among others). Then I passed the sortmerna'd read through FastP to filter out adapters, paired read correction, and low complexity filtering (alongside whatever else FastP does by default). I had ran FastQC before FastP to check if the rRNA filtering did anything and it hadn't, so I also asked FastP to filter by length (--length_required 151) since there were small amounts of reads with lengths shorter than this value which I thought could be the source of the additional peaks in the data.
However, even after sortmerna and FastP together, nothing has really changed:

(There are only two traces/curves here instead of four because the paired reads were apparently recognized properly by FastQC now.)
My plan was to run Khmer to clean the reads a little bit more and then assemble them. I did do a quick assembly with Trinity (with the raw reads) and the BUSCO scores were pretty bad (only ~ 25% of the arthropod BUSCO genes present in the assemblies). I was hoping I'd be able to get Dammit to annotate whatever was left after the cleaning, despite the poor BUSCO prognosis.
Now I'm not entirely sure what to do if I can't clean the reads to begin with, since preventing nonsensical contigs would be absolutely important for a decent annotation.
So this leads me to my questions:

Are these data even salvageable?
Is there anything else I can do to clean the data?
Is this even contamination to begin with that I'm seeing here? (I am presuming some kind of contamination because both samples share that peak at ca. 18% GC content, among other shared features.) Would mapping all the reads against each other and discarding the shared sequences help? If yes, is there a tool I can use for this?
Any other thoughts in this regard?

I searched through the posts here and on Biostars as well as I could but I wasn't really able to find anything helpful in this regard (hence my own post).
I'm very new to all of this, so any inputs, tips, advice, and/or guidance would be most appreciated.

Please can you provide the other metrics given by fastQC? eg the quality per base? I didn't use the latest version by i am suprised by the y-axis in %. In addition, how did you run Trinity? on a merged set of reads or did you perform 1 assembly for 1 sample?

I am sorry but the links to fastqcfiles seems broken. Can you just run a fastp (without filtering according to length, or set a lower size threshold, eg 36 nt) followed by Trinity? You may also give a try to Bridger or BinPacker which are faster than Trinity without quality loss.

also, can you please tell me how you got the data. An academic platform? a private company?

according to your Fastqc files, your reads have a moderate quality. But I am sure there is something to do to improve them.
Try to run fastp with the following options:
fastp -i _1.fq -I _2.fq --detect_adapter_for_pe -l 36 -5 -3 -q 20

and run FastQc again to check to efficacity of the trimming.

sorry I can't see them on the repo... but can you provide the json or html output of fastp too? thanks!

ok thanks. Indeed, no strong improvement. I think you may provide custom fasta sequences to fastp using the --adapter_fasta argument. You may add sequences from Clontech Universal Primer and Clontech SMART CDS kits because contamination most likely come from that. Besides, do you have an idea about the number of genes in the genomes you are interested in? I saw you had less than 10M of reads and i wonder if it is sufficient for a complete de novo transcriptome assembly.

yes you're definitively! Concerning the overrepresented sequences, I never faced this problem for now. But the % is quite high (more than 3%) so i guess it deserved to be tested. As an alternative, you may try to map your RNA-seq reads on the related genome species (use HiSAT2+Stringtie pipeline). It might help to clean your reads!

for sure, do not hesitate to contact me again. Hope you will find a solution for your data anyway.


text:
What are the state-of-the-art cell-type RNA-Seq deconvolution methods?
I would like to find the proportion of each cell-type in bulk RNA-Seq transcriptomics data. I am looking for some guidance on the following:

What are the state-of-the-art methods?
What are their assumptions? (e.g. do they assume that cell-type expression signatures are available in advance?)

Can you describe your experimental design further? This is often solved with single cell RNA-seq

It's not my experiment - I am just trying to deconvolve bulk RNA-Seq data from a publicly available database (GTEx), e.g. given a whole blood sample, find out the proportions of each cell type (e.g. lymphocytes, macrophages,...). Any suggestion?

I see a lot of people using xCell. There should also be papers systematically compare many different methods, which might interest you.

Thanks for the pointer. Do you have any reference of a good up-to-date review paper? A simple description of the methods and requirements would be highly appreciated too

This one looks promising although it may not cover all the methods https://www.biorxiv.org/content/10.1101/2020.01.10.897116v1.full.pdf


text:
Biological datasets for polythetic classification
Question: What are some instances of polythetic datasets in biology? In particular, I am looking for a dataset to benchmark a machine learning algorithm optimized via episodic training.

On polythetic classification
The conventional definition of a conceptual or monothetic class is that its members must possess certain properties in common [1]. For example, in the MNIST database of handwritten digits, images can be classified into 10 distinct classes (i.e. digits) where all instances of any given class possess the same set of attributes (e.g. zeros are typically written as ellipses).
In contrast, a polythetic class encompasses members that might share a number of characteristics, but none of these attributes is essential for membership of that group or class. For example, one could group green zeros together with red ones in one polythetic class, and red zeros together with green ones in another one, as shown below.

Are there any instances of biological datasets with similar properties? I am particularly thinking of genetics data (e.g. single nucleotide polymorphisms) where perhaps a certain mutation A counteracts the effect of mutation B? Any pointers are welcome.

Hi @rvinas, this question seems far too general to be useful. There are, for example, plenty of examples of gene sets that have negatively correlated expression under different conditions, and heritable syndromes that have different physiological appearances depending on the parent they were inherited from. Could you please be more specific about what sort of data you're looking for.

I may not be understanding exactly what you are talking about, but it seems to me that nearly any set of classes of interest must be polythetic. Most groups of comparable things differ in some characteristics but also have characteristics in common, and none of those characters is 100% identical with membership in the class.
Some examples off the top of my head:

Biological species / taxonomic groups (classify according to phenotypic characters)
Biological species / taxonomic groups (classify according to molecular sequence e.g. DNA)
Cell types (classify according to cell epitopes)
Cell types (classify according to transcriptional profile)
Ecosystems (classify according to species profile)
Disease state (classify according to any one of many molecular profiles)
Members of population groups (classify according to genetic variation)

etc.


text:
Challenging benchmarks for supervised learning on sparse scRNA-seq data
One challenging aspect of modeling scRNA-seq data is data sparsity, that is, scRNA-seq measurements typically suffer from large fractions of observed zeros (i.e. dropouts), where a given gene in a given cell has no unique molecular identifiers or reads mapping to it [1].
What would be an appropriate dataset to benchmark supervised learning algorithms that operate on single-cell RNA-seq data? In particular, I am looking for challenging datasets (preferably with labels, i.e. cell type classification) where existing methods struggle due to the sparsity of the data. Additionally, does it ever happen that the dropouts are confounders of the labels? (i.e. for a given cell type A the dropout probability of a gene X is higher than for another cell type B).
Any pointers are welcome.

Can you say what is the purpose of your supervised learning method? to annotate cell types?

Yes, that's right, cell type annotation is the first thing that comes to mind, but I am also willing to consider other tasks (e.g. classification of cancer cells). Even regression tasks would be fine too. Essentially, I have a method in mind that might be useful for modeling non-linear+sparse data and I am looking for an appropriate benchmark dataset.

You could try the PBMC 3k data from the Satija lab:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
For well-annotated data, there's the Single-cell proteo-genomic reference map:
https://www.biorxiv.org/content/10.1101/2021.03.18.435922v1
Data for that are available here.
If these datasets are too large (i.e. too good), you can use the existing counts as weighting for a random removal of reads. Genes with high expression are likely to be present regardless of sequencing depth, hence the count-adjusted removal.
FWIW, I don't think that supervised learning is a good fit for dealing with data sparsity in cDNA read counts. There are plenty of existing attempts to solve the issues by adjusting the assumed models, including research by the authors of commonly-used bulk differential expression tools (e.g. here).

Thanks! I like the second reference, it is fresh from the oven. If the data is "too good" (i.e. traditional models can easily distinguish the different cell types), you would then introduce random dropouts to increase the sparsity?

Random dropouts at the read level; I wouldn't randomly drop out genes.


text:
Do you know a program/script to count how many sequences have mismatches on each position of my primers?
I'm working with primer analysis and design, and so far I've performing it in Excel but: a) I don't like Excel, b) I do it manually, c) I don't like Excel, d) it gets tedious when working with thousands of sequences and tens of primers. My analysis consists of establishing how many sequences have mismatches on each position of my primers, and also count how many mismatches each sequence has with the corresponding primers (see attached image). This last analysis I've been able to (partially) solve it by using ncbi-blast+ tool. However, I will find very helpful a program/script where you input your sequences and primers fasta files and retrieve the number of mismatches at every primer position. Does anyone know something that could help me achieve all this? If it also has an output tab-separated file as the one attached that would be perfect for easy visualization of results. Suggestion/orientations are also grateful! I just want to know if there's something to help me out, otherwise I'll try writing a script to automatize my analysis. Thank you!

Its important to note this isnt a very clear question, for example how many loci are you examining, <10 or 100s? We tend to like knowing the organism and goal of the research. Anyway, at a guess I would just use NCBI's primer-blast program, here https://www.ncbi.nlm.nih.gov/tools/primer-blast/ . This appears to do everything on your list of requests

Michael, thanks for your comments! At the time, I'm working with >50000 sequences and 5 pair of primers/probes for one gene. Primer-blast does not support IUPAC nucleotide nomenclature, so that's a first limitation I've found. Also, I don't find the output information from primer-blast very helpful since it does not summarize the number of mismatches per primer position, and that's what I need in order to improve primer/probe design (I need a count of sequences rather than a list of them). Or perhaps I'm not setting correctly the output format parameters. Any further suggestions? Bests!

I'm amazed at your tenacity in getting Excel to do what you want; that looks like many hours of work to count mismatches.
The hammer I use at the moment for this sort of stuff is LAST, which will do local alignments of sequences to references. If the number of mismatches are not too large (as would be expected for primer sequences), a local alignment will extend to the end of one of the sequences, and I expect that should work fine for your purposes.
I'm going to demonstrate two options, one that can be done just using Excel and web services, and one that needs command-line access. Note that these are only two ways to fetch mismatch counts; you might want to do more than that (it's not clear from your question), and there are many other ways. Also, this may not work for your particular circumstances, but I think it should be quick enough that it's worth a try:
Option 1 - Web-based LAST
In this, I use the Web version of LAST to generate alignments for the primers [as "reference"] mapped to the sequences [as "query"].

To help create the mismatch (or in this case more correctly, "edit distance") statistic, in the Advanced Parameters section, set the a: gap existence cost to 1, and the b: gap extension cost to 1.
Click on Align!
Once the Alignment Results appear, click on the Alignments in tabular format link down the bottom, and save the file to your computer as last.tab.csv (the .csv at the end is the most important bit, to make sure that Excel recognises it as something it can open).
Open that file in Excel, set the delimiter to <tab>
The first column should be the alignment score, which in this case is the match score. This can be converted into a mismatch score by subtracting it from the primer sequence length (column 6).

Option 2 - Local LAST Installation
Installation
Install LAST and samtools.
Index Generation
The first step is to generate a reference index. LAST doesn't care too much about what's query and what's reference, but because I want all the genome sequences to appear on each line, the genome sequences will be the query, and the primer sequences will be the target. I use Martin Frith's lastdb mantra to do this:
$ lastdb -uNEAR -R01 primers.fa primers.fa

Mapping
The next step is to do the actual alignment with LAST. I set the match and mismatch to be both 1 (i.e. the default for the Web version), set the gap existence and extension to be 1, change the output format to tab, and save the result to last.tab.csv:
$ lastal -r 1 -q 1 -a 1 -b 1 primers.fa sequences.fa -f tab > last.tab.csv

Aggregating data
Finally, it's a matter of aggregating the count information. I use grep -v to exclude lines that begin with a #, use awk to subtract field 1 from field 6, then sort the results and use uniq -c to count the mismatches:
$ grep -v '^#' last.tab.csv | awk '{print ($6 - $1)}' | sort | uniq -c
2118 0
   5 1
   5 2


text:
Modify how a script prints the output?
I'm trying to print out all measured distances between a ligand and 2 residues (using PyMol), below is the script that I'm using. I only want to measure the distances between all the ligand's O & ND2 of N312, as well as all the ligand's N & OD2 of D113.
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "/r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("dst4.txt", 'w')

for a in stored.a:
  for b in stored.b:
    for c in stored.c:
      for d in stored.d:
        f.write( "distance from %s to %s is %s\n" % (a, c,
cmd.distance( "id %s" % a, "id %s" % c,cutoff=5,mode=2)))
        f.write( "distance from %s to %s is %s\n" % (b, d,
cmd.distance( "id %s" % b, "id %s" % d,cutoff=5,mode=2)))

f.close()

The output file currently looks like this
distance from 9 to 1334 is 4.589359760284424
distance from 8 to 4204 is 3.1509475708007812
distance from 9 to 1334 is 4.589359760284424
distance from 15 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 17 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 20 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 8 to 4204 is 3.1509475708007812
distance from 23 to 1334 is 0.0
distance from 15 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 17 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 20 to 4204 is 0.0

I want my output to look more like this. Thanks in advance!
Lig_N D113 4.59
Lig_O N312 3.15

Let's simplify...
cmd.distance(xx,xxx) accepts two selector strings and is fine with and state 1. And best of all if you select multiple atoms it gives you the smallest distance!
So all you need to do is:
# measure
d_O = cmd.distance('UNL and element O', 'resi 312 and name ND2')
d_N = cmd.distance('UNL and element N', 'resi 113 and name OD2')
# write
f=open("dst4.txt", 'w')
f.write(f'Lig_N D113 {d_N:.2f}\n')
f.write(f'Lig_O N312 {d_O:.2f}\n')
f.close()

If you are using python 2 (if so do note pymol scripting is better in 3), then the second bit is
f=open("dst4.txt", 'w')
f.write('Lig_N D113 %.2f\n' % (d_N,))
f.write(f'Lig_O N312 %.2f\n' % (d_O,))
f.close()

Unsolicited coding tips
A few unsolicited tips... Your original script is rather over-complicated and obfuscated and uncommented, so for better legibility why not do:

cmd.iterate_state is always harder to follow than the for atom in get_model('xxxx').atoms synthax
unless you are stuck with python 2, use f-strings not the %s sprintf style formatting
the /x/x/x/x selector synthax is never as readable as the longer
resi/resn/name selection algebra.

Code golf is always tempting, but it can result in future-you getting stuck as happened here...

You're right about the unnecessary complication in my script haha... I adapted the script a bit and it works :D Thank you!


text:
RNA-Seq data transformation prior to sample correlation analysis
If I’m starting with Deseq2 normalized counts what are some preprocessing steps that I should apply to these data before estimating sample correlation using the cor function in R? For example, would it make sense to quantile normalize the data?

Hi @Alexis your question doesn't quite make sense, you are starting with normalised counts and you want to further normalise them? If you data is normalised you can proceed with Spearman's correlation analysis.

You don't absolutely need to perform any further steps before calculating the correlation. However you would likely benefit from the following sort of workflow:

Transform your data with vst() or robust log (see the DESeq2 vignette) so genes are on a more or less similar scale.
Select the top 300-500 most variable genes
Use this subset with to compute sample-wise correlation.

This also what follows what you would find in the DESeq2 user guide, which I strongly encourage you to read.


text:
MM/PB(GB)SA calculations for DESMOND trajectories
I did molecular dynamic simulations with DESMOND.
Now, I want to obtain the interaction free energy using MM/PB(GB)SA.
The free tools I found to perform such calculations are AmberTools and g_mmpbsa.
AmberTools work only on Amber-type files, whereas g_mmbbsa works with GROMACS-type. This imply that I have to convert the DESMOND trajectories to AMBER or GROMACS formats.
I did that using VMD. In case of AMBER, I got the files prmtop and mdcrd whereas for GROMACS I got gro and top files. I also tried to use InterMol, but I am getting an error on conversions.
For both tools, I need to prepare the input files. In case of AmberTools, there is a Python script (ante-MMPBSA) that extract the solvated complex, the complex, the ligand and the receptor from the trajectory files. This scrip is falling to do that from the converted trajectories. In case of g_mmpbsa, it also need the index.ndx file.
My question: Is there a way to accomplish this task? Any suggestion is welcome!

Try here: https://mattermodeling.stackexchange.com/. There's way more VMD questions there, plus there's an entire tag just for GROMACS, and AMBER is on topic there too.

Hi @user1271772 (or I must said @Nick? This is Camps, lol.


text:
How to get a file with the introns of a gene?
I was able to get the exons of a particular gene as a fasta file from the NCBI genome viewer. As for introns I only found a bed file that showed where the intron starts and ends, along with its direction. What tool do I use to get the exact sequence of these 11 introns using this bed file?
And is there a more direct way?
Thank you.

First I edited the Chromosome column of the bed file from ncbi, because it contained the accession number instead of the chromosome. Then I used the getfasta tool of the bedtools utility.
bedtools getfasta -fi input_genome_file.fa -fo output.fa -bed bed_file
The output contained the introns in the format I wanted.
If anybody had a better way, feel free to post.

Hi @venkateshwar thank you for sharing this information. Could you just confirm that this is the answer to the question you posted (it looks like it)?

@M__ Ok I will.


text:
How to reverse complement a genbank file?
I am visualizing the alignment of a particular genomic region across different genomes with mauve. In certain genomes the region of interest is in the negative strand and in others the positive strand. What tool can I use to reverse complement my genbank files without losing the features and annotation?

I have not really used Mauve, so I don't know whether it might have an in-built setting for that.
SnapGene Viewer should be able to do that - it can open genbank files and export them, and View > Flip Sequence gets you the reverse complement without affecting the annotations. SnapGene can also do that; I only specified Viewer because it's free.

Thank you. It worked.

From the biology side of bioinformatics and the easier option, I would use Mauve.  If you want to make the project more specific and tailored to your needs, and you have python skills, I would write Python code to complete this task.  I would specifically use the bioinformatics Python library called Biopython: https://biopython.org
pip install biopython
The second option is more involved but it gives you more control on your project/research since you're the one creating it.  I personally would go with the 2nd option if you have coding skills.  If you don't have coding skills, stick to the 1st option.  Im sure there are also tools you can use that are found on Galaxy:
https://usegalaxy.org
I am not sure what tools there are that would be good for your task and needs, but from using Galaxy in the past, im sure you can find something! :)


text:
Calculate Entropy for DNA Multiple Sequence Alignment in R
I am pretty new to R. So I apologize for asking maybe a very basic question.
Let's say I have a fasta file with sequence below:
>sequence_1
ACCTGC--A
>sequence_2
ACC-GCTTA
>sequence_3
ACCTGCTTA

Is there a function or method to calculate entropy for each column in R?
I did some research and found a function that calculate entropy for Amino Acids but not for DNA nucleotides.
I would like to have the output be stored in a vector or in a list.
Thanks for the help in advance!

Hi what is the goal of your research? Entropy is a function of a given allele since you have variation, particularly indels this is not possible. With an indel position e.g. 4 now becomes position 3 (sequence_2) within an entropy calculation, so which column are you referirng to? I think you are confusing this with phylogenetic methods.

Hi Michael, I was able to get it to with some other people's help. But to answer your question, The MolecularEntropy function in R did exactly what wanted to. It calculated the entropy at each nucleotide position with aligned sequences.

okay, the blurb does say "Gaps are ignored on a site by site basis so the entropy values may have different number of observations among sites." It also needs sequences to be the same length ... I personally don't think it works with indels


text:
How to retrieve the DNA sequence from a particular species and region through command line?
I am new to working with the NCBI database so I am not sure how trivial this question is. I wanted to get DNA sequence of a specific version of the Xenopus tropicalis genome of a specific region (ie 1-100bp) without doing it manually each time (going on genome browser). I know there is a faster way and efficient way to do this through the command line but I found the NCBI tutorial hard to follow.

Hi @fp1234 perhaps you could post the code you have tried to date? I presume you are using the efetch system? There are other alternatives.

Hi! Yes I am using efetch. So as a test I was trying to run  efetch -db genome -id 5323661 but I keep getting a 400 bad request error. The id is the UID for this genome version Xenopus tropicalis https://www.ncbi.nlm.nih.gov/assembly/GCF_000004195.4/. But would I hope to do is to be able to get "x - y bp" from z chromosome using efetch.

Samtools
Samtools is a toolkit that includes a binary called samtoools, which can be used to query a FASTA file — to extract a subsequence based on some coordinate range you specify.
If you do not have samtools installed, you need to add it to your computer. For OS X, you could use Homebrew, for example:
$ brew install samtools

For Ubuntu, you could use apt-get:
$ sudo apt-get install samtools

Reference assembly
The next step is to download a reference assembly for your genome and genomic build-of-interest.
There is a research group dedicated to making Xenopus datasets available. If you were interested in v10 of the Xentr genome, for instance, you could do the following or similar to download it:
$ curl --output "XENTR_10.0_genome.fasta.gz" "http://ftp.xenbase.org/pub/Genomics/JGI/Xentr10.0/XENTR_10.0_genome.fasta.gz"

You could also dig through the NCBI site to get a similar link. It's the same idea — just change the web address, accordingly.
Indexing
Indexing makes it possible to extract an arbitrary subsequence of interest from this FASTA file. You need to index, but you only need to do it once:
$ samtools faidx XENTR_10.0_genome.fasta.gz

This creates two new files: XENTR_10.0_genome.fasta.gz.fai and XENTR_10.0_genome.fasta.gz.gzi.
Note: You need to keep these two files in the same directory as the original, compressed FASTA file XENTR_10.0_genome.fasta.gz.
Querying
Now you are ready to query, or extract a subsequence of interest. If you wanted the first hundred bases of the first chromosome, for instance:
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100
>Chr1:1-100
cctcagcgtcttcttctttctgctcccactgcattagcatagggagagagggccacagag
gctaggtagcattgccggctagcaatcttcagcgtcctct

The search term used here is Chr1:1-100. Replace Chr1 with the chromosome of interest and 1 and 100 with the start and stop coordinates, respectively, of interest. Coordinates used with samtools queries use a one-based indexing scheme.
Note: This uses the chromosome naming scheme in the XENTR_10.0_genome.fasta.gz file (Chr1, etc.). Other research groups (NCBI) might use different chromosome names. If you get your FASTA from somewhere else, investigate the FASTA file to determine the naming scheme being used.
If you want to redirect output to a file, use the output redirect operator (>):
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100 > mySubsequence.fa

I'd download the fasta, and use samtools, or maybe bedtools.


text:
How can I restrict the search space when querying GEO?
I am using the following command to retrieve GEO Datasets entries for daf-2 c. elegans that have as their experiment type High Throughpout sequencing:
esearch -db gds -query "daf-2" | efilter -query "expression profiling by high throughput sequencing [DataSet Type] AND Caenorhabditis elegans [organism]"

However, among the results that I got, I noticed that I got this:
23. Impaired Insulin-/IGF1-Signaling Extends Life Span by Promoting Mitochondrial L-Proline Catabolism to Induce a Transient ROS-Signal
(Submitter supplied) Transcriptome profiling of three models with impaired insulin/IGF1 signaling. 1. Deep sequencing of endogenous mRNA from Caenorhabditis elegans N2 var. Bristol (wildtype) and daf-2(e1370) mutant; 2. Deep sequencing  of endogenous mRNA from murine embryonic fibroblasts (MEF)  wildtype and irs1-/- knockout; 3. Deep sequencing of endogenous mRNA from murine embryoinic fibroblast (MEF) insr+/- -lox and insr+/- knockout  Jena Centre for Systems Biology of Ageing - JenAge (www.jenage.de)
Organism:   Mus musculus; Caenorhabditis elegans
Type:       Expression profiling by high throughput sequencing
Platforms: GPL11002 GPL13776 14 Samples
FTP download: GEO (CSV) ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE36nnn/GSE36041/
SRA Run Selector: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA151765
Series      Accession: GSE36041 ID: 200036041

As you can see, in the Organism field, there are both Mus musculus and Caenorhabditis elegans. How can I restrict my search so that I only get entries that exclusively work on Caenorhabditis elegans? One naive way to do this would be to include NOT Mus Musculus in the query, but of course that would mean that I could get something else too.
Another way I can think of solving this is just to write a script that would do the extra filtering using regexes, but I was wondering if there's a simpler solution using the e-utilities functionalities

That particular study was done using both worm and mouse samples. If you are ultimately interested in SRA data, you should just pipe the output of this step to elink and get SRA run table. You can either add an organism filter after that or just subset the run table.

@vkkodali but why would filtering from SRA's side be different from filtering from GEO's side?

@vkkodali I noticed that it does filter out other organisms, but I just don't understand why filtering from SRA is different from GEO. Also, I'm filtering again by "Expression profiling by high throughput sequencing [DataSet Type]" because of the same issue as with organisms, but when I do this, I get 0 results. Any clue as to why this might be happening?

Filtering from SRA would filter it at the run level where as filtering in GEO at the series level will find any series that contain at least one sample that is from worm.

@vkkodali I see that makes sense. What about the other filter?

You can use the following query to first search GEO to get a list of GEO Series of interest, then find all linked SRA runs and perform a second filter step to keep only the SRA runs that satisfy another set of criteria as shown below.
esearch \
  -db gds \
  -query "daf-2[all fields] \
    AND expression profiling by high throughput sequencing [DataSet Type] \
    AND Caenorhabditis elegans [organism]" \
  | elink \
  -db gds -target sra \
  | esearch \
  -query "(#3) \
    AND Caenorhabditis elegans [organism]" \
  | efetch -format runinfo

Here, I am combining two independent queries using the (#3) term in the second query. More information about this is in the "Combining Independent Queries" subsection of "Searching and Filtering" section of the Entrez Direct documentation.

could you explain why you have a -db gds with elink?

It’s just a habit. If you are piping results from another edirect tool, you don’t need to specify the db. But if you are running elink as a stand-alone command then you need the db.


text:
Can I use a tool for highly variable gene selection for selecting highly variable "bins" from an ATACSeq dataset?
I'm working with ATACSeq data from multiple tissue/cell types. The data is binned in 1 megabase bins. I'd like to identify the bins that are "highly variable" across the different tissue types. I can't seem to find a tool to do this, like there exists for identifying highly variables genes in a gene expression setting (e.g Scanpy's scanpy.pp.highly_variable_genes). I was wondering if I could use the gene selection tools for my data, and if there are any other suggestions?

I don't think is a good idea to apply the gene selection method on the ATAC-seq genomic regions because of the high sparsity of ATAC-seq data.
A common way is to select the most frequent regions (like the Signac does), or to discard the house-keeping regions (like the SnapATAC does).
And you can also use the episcanpy, a epigenomics version of scanpy.

Thanks for your answer. I guess in my case the cells are grouped per cell type (the data is aggregated per cell type), so sparsity shouldn't be an issue. Given this, is it still inappropriate to use gene selection methods?

If you think the sparsity is not an issue for your data, then scanpy.pp.highly_variable_genes may be applicable. Just be aware that the default thresholds for mean and disp are for single-cell RNA-seq data.


text:
Programmatically retrieve Metadata from SRA Run Selector
I previously asked a question about how to retrieve the Accession List associated with a SRA project. The answer was:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo
where PRJNA491191 is the bioproject that I'm interested in in this example. How can I do the same for the metadata file (the arrow in the picture below is pointing to accession list. It should point to metadata, which is what I'm interested in now)?

Depending on exactly which information you seek, you can try:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format native -mode xml

This will return an XML with additional information that is not present in the runinfo table. You can use the Entrez Direct tool xtract to parse the output XML and extract the fields of interest into a tabular format.
Alternatively, you can try:
esearch -db sra -query 'PRJNA491191[bioproject]' | esummary

to retrieve a summary output that includes fewer fields than the efetch option shown above. Perhaps all of the fields that are of interest to you are present in the summary output... it will save you time and effort in parsing the XML and retrieving the information you want.
For example, you can use the following to retrieve some fields using xtract from the esummary output:
$ esearch -db sra -query 'PRJNA491191[bioproject]' \
  | esummary \
  | xtract -pattern DocumentSummary -element Summary/Title,Statistics@total_size,Study@acc
GSM3391430: hsf-1(OE) rep3; Caenorhabditis elegans; RNA-Seq  731402855  SRP161795
GSM3391429: hsf-1(OE) rep2; Caenorhabditis elegans; RNA-Seq  703822418  SRP161795
GSM3391428: hsf-1(OE) rep1; Caenorhabditis elegans; RNA-Seq  715482490  SRP161795
...

Thank you. So I want to create a csv file from the output. Also, I'd like to keep all columns. How can I do this? I'm trying esearch -db sra -query 'PRJNA491191[bioproject]' | esummary | xtract -pattern DocumentSummary -sep "," -element but I'm not sure what should go after -element since I want to keep everything

Also, I can't specify all columns because I don't know what they are beforehand (I know I can look them up from the output but I'm doing this for many projects so I practically can't do this manually)

You have to specify the elements to be able to retrieve the fields you are interested in. Unfortunately, you cannot retrieve everything using xtract without knowing the field names. I will update my answer above with an example.


text:
getGeneLengthAndGCContent in EDASeq complains about "Ensemble gene ID"
I am trying to use the getGeneLengthAndGCContent function from the EDASeq library to retrieve gene lengths for c. elegans. I retrieved Ensemble gene ids from Biomart for c.elegans and I'm using those as input to getGeneLengthAndGCContent (I picked "ensemble_gene_id" from the filters of the worm annotation Mart). Here's an example of how I'm trying to use the function:
getGeneLengthAndGCContent("WBGene00001042", org, mode=c("biomart", "org.db"))

When I try to run this, I get
Error in getGeneLengthAndGCContent("WBGene00001042", org, mode = c("biomart",  :
  Only ENTREZ or ENSEMBL gene IDs are supported.

@ATpoint oh that's weird why does BioMart list it as such?

@ATpoint also when I look at Ensembl I don't find any other identifier for a gene of c elegans (here's an example http://uswest.ensembl.org/Caenorhabditis_elegans/Gene/Summary?g=WBGene00001042;r=III:5415556-5419565;t=W03A5.7.1). Could you guide me as to where to find those other ids?

@ATpoint so WBGene00001042 is indeed the Ensembl id?

I found my answer here. It uses getGeneLengthAndGCContent but changes a small part of the function.


text:
What is a good RNA seq normalization method that allows for across sample comparisons and between transcripts
What is a good RNA seq normalization method that allows for across sample comparisons, and allows between transcripts comparisons as well? I read that TMM for example allows across sample comparisons but not between transcripts. On the other hand, I read that TPM is not good enough for between sample comparisons

I found a recent technique that allows for both, called geTMM.

I'm having a bit of trouble understanding the question, but we've used DESeq2 for doing both between-sample and between-transcript comparisons.
When comparing transcripts for short reads I use the variant-stabilised transformation, adjusted for transcript length. This is good for a visual overview of results (e.g. in a heatmap).
When comparing samples I use the variance-adjusted DESeq2 test, producing normalised $\textrm{log}_2FC$ values. This works well for quantitative results.
See this paper, where we use both of these approaches (e.g. Figure 2 shows both in a side-by-side view).

As mentioned by @gringer DESeq2 is a popular choice for normalization of counts between samples. However, because each gene is adjusted individually it is not well suited to compare within a sample.
For this a global normalization like TMM from edgeR would be more appropriate.
You can find more information in this tutorial from the hbc at Harvard.


text:
What does the "filters" argument in "getBM" do?
In the biomaRt package, there is a function getBM which among a few things is useful for mapping between different gene representations. Currently, I'm using it as follows, to map Ensembl transcript IDs to gene names and Ensembl gene IDs.
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       values = as.list(transcripts),
                       mart = mart)

How is the following code different:
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       filters = "ensembl_transcript_id",
                       values = as.list(transcripts),
                       mart = mart)

It's unclear to me from the description of the filters argument here. But just based on results, I get less results when I run the second one.

Hi An Ignorant Wanderer,
Filters allow you to narrow down the dataset to query a specific subset of features within your dataset. A description and tutorial of the filtering steps within BioMart can be found here:
https://www.ensembl.org/info/data/biomart/how_to_use_biomart.html
If you remove filters, you will query the full dataset (all human genes, in the example above). However, this is not recommended for BioMart as large, genome-wide queries will overload the BioMart servers. For queries of this size, we suggest using the Ensembl REST API or FTP site.


text:
Finding discriminating positions for a sequence alignment column
I don't have any theoretical background in biology, so forgive me if my question is a bit.. well..dumb. I'm trying to use a Monte Carlo approach to find the discriminating position of a given sequence alignment column. So essentially, given a sequence alignment, I want to find the position at which we can distinguish active and inactive conformations.
So, here is my code for the Monte Carlo method:
def MonteCarlo(S, max_iter, threshold):
    """
    Computes the discriminating position for a sequence column

    Input: S = sequence column, max_iter = max iterations, threshold = threshold of acceptance for boltzmann distribution (between 0 and 1)
    Output: 2 classes
    """
    n = len(S)
    S_copy = S.copy()
    num_class1 = int(0.30*n) #number in class1
    print(num_class1)
    for _ in range(num_class1): #intialisation
        rand_position = random.randint(0,len(S)-1)
        if S_copy[rand_position] == 'X':
            rand_position = random.randint(0,len(S)) #reselect if used position
        S_copy[rand_position] = 'X'

    #Monte Carlo with Insertions and Deletions

    s_entropy = entropy(S)
    n_iter = 0 #counter for iterations
    class_1 = dict() #class 1 determined by X's
    class_2 = dict() #class 2 is rest of S
    for i in range(len(S)):
        if S_copy[i] == 'X':
            class_1[i] =S[i]
        else:
            class_2[i] =S[i]

    res = score(S, class_1, class_2)

    while n_iter < max_iter:
        r = random.random()

        if r < 0.25: #25% chance to do an insertion
          (c1,c2,pos) = Insertion(S, class_1, class_2)

        elif r < 0.5: #25% chance to do a deletion
          (c1,c2,pos) = Deletion(S, class_1, class_2)

        else: # 50% chance to do a swap
          (c1,c2,pos) = Swap(S, class_1, class_2)

        res2 = score(S, c1, c2)

        #Boltzmann Distribution coefficient
        delta_S = abs(res2 - res)
        lambda_, N = 0.2, 19
        size = 6000
        ind1 = random.randint(0,size - 1)
        rv = boltzmann.rvs(lambda_, N, size=size)/N
        beta = rv[ind1] #pick float between 0 and 1 with boltzmann distribution
        if  (res2 > res) or (beta*delta_S > threshold): #if better score i.e. higher entropy gain
            res = res2
            class_1 = c1
            class_2 = c2

        n_iter += 1

    return class_1, class_2, pos

Here I'm using the Boltzmann distribution as a probability measure. All the other functions (swap, insertion, deletion, entropy) calculate pretty much what the name says (if you need more clarification on my side, I'd be happy to provide it!)
Let me explain what my code is trying to do: Given a sequence alignment (as a list) say ['M','M','M',...,'L','L','L'] I first set a bunch of random positions and mark them out as 'X' (the number of Xs are determined by num_class1). Next, in the while loop I do swaps,insertions, and deletions with the goal of obtaining a better score. This meansI'm adding, removing, or swapping positions between class1 and class2. Eventually, my code returns two classes and the discriminating position.
While my code works, I don't think it's giving the right result. When I use a base case example of a list with two letters 'M' and 'L with a clear discriminating position, I get random results.
Could someone explain theoretically (or via code) how one can find the discriminating position using a Monte Carlo method as above?
Thanks!

It's not very clear what you mean by "discriminating position". Can you explain more or link to a source that describes in a little more detail what you are talking about?

Hi! By discriminating position, I mean I'm trying to find a position within the sequence alignment that splits the columns into two classes such that the two have the lowest entropy possible. So, in my sample case with a clear division between Ms and Ls, I should've gotten the the position at which the Ls begin. Does that make sense?
Here is a link: https://pubs.rsc.org/en/content/articlelanding/2019/sc/c8sc03669a#!divAbstract

Hi, the paper that you've linked to is talking about a discrimination between two structural states of a protein by a kinase; "discriminating position" does not appear to occur in its full text. It's not clear how this generalizes to alignments. From your explanation though I think I have a better sense of what you mean: for seq alignment $S$ with positions i = 1,2,3...n find the position $S_i$ such that a measure of entropy across aligned positions $S_1 : S_{i-1}$ and $S_{i+1} : S_n$ is minimized for both slices across each column of the aligned slice. (words are easier than code for me).

Yes! no problem (: I also misunderstood what I was trying to solve actually, I have managed to fix the problem! But for future reference: the code is supposed to split the sequence column into two classes such that each distinctly has the highest entropy possible (or lowest for that matter)

Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.

Right I see, so I need to preprocess the sequence alignments first? So that I match the gaps properly?

Yes, to verify this step for logical errors


text:
How do I pull singe cell RNA sequencing data from GEO database?
I am new to R and computational biology. I am trying to look through a published data set to check for gene expression for my own project. I am having trouble finding materials to teach me how to analyze single cell RNA sequencing data on R from the GEO NCBI database.
This is the data that I am trying to access https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138651
I would appreciate any directions or suggestions for tutorials.
Thanks!

Install package
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GEOquery")

Use function getGEO to download the files.
expressionSet <- getGEO("GSE138651")
expressionSet <- expressionSet[[1]]

You can then look at the info of the samples:
pData(expressionSet)

And the assay data itself
exprs(expressionSet)

This is useful because of the very specific instructions but I think it would be even more helpful to also include a brief introductory explanation and links to BioConductor and the GEOquery package to provide context for newcomers.

I am sorry but I have to flag this for removal since it's wrong. getGEO only returns expression data for arrays, not for NGS. Output of exprs(expressionSet) is empty and this is expected.

The simplest would be using a count matrix (at the end of the link you have shared, section "Supplementary file"). For example GSE138651_Vagal_WholeNodoseSeq_raw.mtx.gz corresponds to "raw" count matrices obtained after mapping and counting the reads. This file, along with the GSE138651_barcodes.tsv.gz and GSE138651_genes.tsv.gz can be used to "read" the data into R; to that end there are many options in the form of packages, Seurat is one of the best documented and hence is beginner-friendly.

Thanks Haci. I am learning Seurat as well. So many things!!
I am actually the most interested in the file GSE138651_NGTall_Reads_cap.csv.gz. Do you know if this file type can be "read" into R?

Yes, it is a csv file and once again there are many options to read such files into R. read.csv() is one option. Personally, I use fread() of data.table to read in tabular files.

I will try that then. Thank you again!

You might want to check out the BioConductor GEOquery package.
Install using
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GEOquery")

More information on how to use this package and read in data from GEO can be found in the vignette

I am sorry but I have to flag this for removal since it's wrong. getGEO only returns expression data for arrays, not for NGS. This package does not apply here.


text:
How do I generate a list of post-synaptic gene markers that will guide me in my search for these markers in intestinal stem cells?
I would like to generate a list of gene markers that represent post-synaptic related genes. I will then use this list to search through single cell sequencing data, in order to assess whether certain cell types express post-synaptic related genes. I think there might be a way to get this list of gene markers from NCBI.
Thank you!

something like this ? https://www.gsea-msigdb.org/gsea/msigdb/geneset_page.jsp?geneSetName=REACTOME_NEUROTRANSMITTER_RECEPTOR_BINDING_AND_DOWNSTREAM_TRANSMISSION_IN_THE_POSTSYNAPTIC_CELL

Yes! i think so. Thanks for sharing. I didn't know about this database!


text:
Output when compute GC content is different between using user input and using string
I am learning bioinformatics and i was solving a problem on rosalind.info.
Input is a string consist of multiple DNA strand names and the DNA sequences, and Output is the name of the sequence that has the highest GC percentage and the GC percentage.
this is the input:
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

this is the correct output:
Rosalind_0808
60.919540

This is my code (sorry for this filthy mess):
DNA = input()
wholeRosalind_DNA = DNA.splitlines()
Rosalind = []
Rosalind_DNA = []
for i in wholeRosalind_DNA:
    if "Rosalind" in i:
        Rosalind.append(i)
Countstring = []
maxGCcontent = []
GCcontent = []
def eachRosalind_DNA(DNA):
    for i in range(0,len(Rosalind)):
        if Rosalind[i] in wholeRosalind_DNA:
            a = wholeRosalind_DNA.index(Rosalind[i])
            if a == wholeRosalind_DNA.index(Rosalind[-1]):
                b = len(wholeRosalind_DNA)+1
            else:
                b = wholeRosalind_DNA.index(Rosalind[i+1])
            Rosalind_DNA.append(wholeRosalind_DNA[a:b])
    for i in Rosalind_DNA:
        Countstring.append(''.join(i))
    for i in Countstring:
        bases = len(i)-14
        CG = i.count("C") + i.count("G")
        GCcontent.append([i[1:14],(100*(CG/bases))])
    for i in GCcontent:
        maxGCcontent.append(i[1])
    result = max(maxGCcontent)
    for i in GCcontent:
        if i[1] == result:
            return (i[0]),(round(i[1],6))
print(eachRosalind_DNA(DNA)[0])
print(eachRosalind_DNA(DNA)[1])

but i get this instead:
  File "GC_Content.py", line 26, in eachRosalind_DNA
    GCcontent.append([i[1:14],(100*(CG/bases))])
ZeroDivisionError: division by zero

i will get the desired output if i just paste the input into the variabe "DNA" without using user input. Please help me fix the code and distinguish those two ways of getting input

Piece of advice: use meaningful names for your variables (loop variables do not need to be named i). This will help you and others spotting potential bugs.

Is it cheating to use Biopython(?), will a solution might be helpful anyway. This assumes you have your sequences stored in the seq.fa FASTA file:
from Bio import SeqIO
from Bio.SeqUtils import GC

print(max(((r.id, GC(r.seq)) for r in SeqIO.parse("seq.fa", "fasta")),
          key=lambda x: x[1]))

Output:
('Rosalind_0808', 60.91954022988506)

I think the hard thing in this problem is not the calculation of GC content, it is the processing of DNA strings because they are separated on different lines. I would suggest reading the file line by line to generate complete DNA strings for each Rosalind entry. (I am not able to spot what is causing the error in the code).
seqDic = {}
with open("YouInFile", "r") as inFile:
    for ln in inFile:
        ln = ln.rstrip()
        if ln.startswith(">"):
            ln.name = ln[1:]
            seq = inFile.readline().rstrip() + inFile.readline().rstrip()
            seqDic[ln.name] = seq
        else
            continue

Basically, whenever you encounter a ">" in your file, you write down the Rosalind name, and then you read the next two lines and combine them together. You put the name and the sequence in a dictionary, and iterate over the dictionary to calculate GC content.
Note that I did not debug my code. Let me know if it is not working.

The error is that bases = 0. Briefly looking at the code carriage returns in your original code will trigger this, as will confusing a faster header with sequence. The latter is because the fasta headers are 14 characters long and subtracting 14 from bases which is of course zero. The former will likely result in 'nan' because it would be 14 subtracted from ''.
I've a simpler solution, download Biopython (using conda or last resort pip), you might want to use SeqIO (Biopython), which would be cooler because it allows you to call Rosalind_0808 directly. Then:
from Bio.SeqUtils import GC
seq = 'CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT'
print (GC(seq))

Done. If you want to splice off the last 14 characters its
GC(seq[:-14])


text:
Using Entrez.efetch() to retrive .fasta file from any NCBI database?
I tried to use these python codes to retrieve a .fasta file from gene database, but the output were the parameters of the report
from Bio import Entrez
Entrez.email = "send_helps@earth.com"
handle = Entrez.efetch(db = "gene",id = "682",rettype="fasta", retmode="xml")
print(Entrez.read(handle))

i also read the documentation on this site about efetch() and seems like this function can only retrieve certain kinds of data for each database. So i can't take fasta file from gene database !?
I don't want to use the NCBI web UI. What's wrong with my code? Is there a way i can guarantee to retreive a .fasta file through python from NCBI ? Am i misundertanding something severly ?
Sorry for bad English. thank you.

In order to download the FASTA sequence of transcript variants using Entrez eutilities, you will have to use Entrez.elink to get a list of transcript accessions for the gene of interest and then use Entrez.efetch to fetch the sequence(s). Just to give you an idea, you can use Entrez Direct for this as follows:
elink -db gene -id 682 -target nuccore -name gene_nuccore_refseqrna \
  | efetch -format fasta

You can find more Entrez Link descriptions here. Note, there is a gene_fasta option for efetch but it still returns only the transcript sequence (i.e., no intron sequence).
If you need the full gene sequence (including the intronic regions), you can use the commands mentioned here.
Alternatively, you can use NCBI Datasets for this. You can use the web interface, command line tool or the python API. This page with the how-to articles has the information you seek.


text:
BigWig circular plot representation
I would like to ask if anyone knows or has experience in representing BigWig signals in a circular plot? I mean like IGV but in circular plot not linear. Or any library that can be used in plotting BigWig signals especially in Javascript?
I appreciate any suggestion,
Thank you.

circosJS should probably work well, based on your requirements: https://github.com/nicgirault/circosJS
I posted a short Observable notebook that shows what calls could be used to bring data from remote bigWig files into a circosJS plot:
https://observablehq.com/@alexpreynolds/rendering-bigwig-file-data-in-a-circular-representation
It is less intended as an explicit how-to in Javascript, than as a starting point for the end user with some of the function calls and data scale issues involved. Hopefully it can be useful to the reader for exploratory purposes.

Does it take bigwig file without any conversion to tabulated format?

Doesn't look like it, but you could probably use https://github.com/GMOD/bbi-js as a file parser, get the data you need asynchronously, and feed it to circosJS to render.

Thank you very much on your suggestion, I will give it a try and see if it works.


text:
Chromosomal order of supplementary alignment in BAM file "SA" tag
I have several long-read sequences, which when aligned to a human reference genome, aligns to multiple chromosomal fragment as a result of chromothripsis shattering the linear DNA and randomly piecing them together. This creates for a single primary alignment multiple supplementary alignments.
With a single example of a read (only showing parts of all the bamfile information) being
cdaf1fcb-27bf-4890-8ac7-a097dfa12c38    16      chr11   15151487        60
With the primary alignment starting at chr11 at 15151487
Then the supplementary alignment can be seen as follows
SA:Z:chr13,49827969,+,1839S956M25D1584S,60,53;chr11,15152459,-,833S749M22D2797S,60,72;chrX,138324656,-,2540S404M2D1435S,60,20;chr11,15152066,+,278M11D4101S,60,20;
Where the alignments are to chr13, 11, X and 11. With both regions on chr11 being close to the initial region from the primary alignment, making the regions on chromosome 11 the same (I know I haven't calculated the full chromosomal interval using the Cigar string for this example).
I assume the order of that supplementary information above is in the actual order. But the problem is that my DNA should be circular which doesn't fit with the order of the alignments since the circle would then be chr 11 (primary, start of circle) -> chr 13 -> chr11 (same as primary) -> chrX -> chr11 (back to the start of circle)
My question is:
Q) So is the order of the supplementary information just pieces together randomly or is it true that I can assume the order is the true order of the sequencing read?
So thanks a lot for your help.

There is nothing in the SAMtags specification that dictates the order of the entries in an SA tag.  It may be entirely dependent on the aligner that you used.
Looking at my own data that was aligned with bwa, I cannot see a definite pattern.  Longer matches tend to be first, but not always.
Your aligner documentation may have more information on how it uses SA tags.


text:
Use Htslib to extract specific region of vcf file
I am working in c++ using the htslib library to parse different vcf files as well as fasta files.
My goal is to read in a specific region from a fastafile and then extract all the VCF sites within that region from a seperate vcf file, in this case only containing chr14. I can extract the fasta region in several ways, one of which is to use the function faidx_fetch_seq() and define a region of chr14:20000000-22000000.
Then for the bcf files i am able to iterate through the entire file and extracting the different information regarding the reference, alternative allele and so on. With a simple example shown below.
void faVCF(const char* fastafile,FILE *fp){
  htsFile *test_bcf = NULL;
  bcf_hdr_t *test_header = NULL;

  test_bcf = bcf_open("chr14.vcf", "r");
  test_header = bcf_hdr_read(test_bcf);
  bcf_hdr_set_samples(test_header, "20000000-22000000", 0);
  bcf1_t *test_record = bcf_init();
  while(bcf_read(test_bcf, test_header, test_record) == 0) {
              std::cout << bcf_hdr_id2name(test_header, test_record->rid) << "\t" <<
                          test_record->pos << "\t" <<
                          test_record->n_allele << "\t" <<
                          std::endl;
          }
}

However, i doesn't work to extract all the variant sites in a specific region. I tried bcf_hdr_set_samples with the argument "20000000-22000000", to extract that given region from my chromosome. As far as I understand this should be able to subsample both samples and region queries using the header information, and then upon initiation extract that region, however this doesnt work.
So my question is: Can anyone help me with what function I can use to extract all the vcf sites from the region chr14:20000000-22000000.?

For iterating over reads in a region in a single BCF file, you want to use the iterator interface: bcf_itr_queryi() or bcf_itr_querys(), together with bcf_itr_next(). For a VCF file, it appears you need to use the similar TBX iterator interface: tbx_itr_queryi() or tbx_itr_querys(), together with tbx_itr_next(). See HTSlib's tabix.c for reasonably straightforward examples of both.
Alternatively you could use the routines in htslib/synced_bcf_reader.h, which encapsulates both and is what bcftools mostly uses.


text:
Snakemake rule issue with bcftools view filtering
I am creating a bioinformatics pipeline using Snakemake doing everying from adapter trimming of raw sequencing data in .fq to mapping, SNP calling and SNP filtering.
The issue i am having is when i try to filter specific chromosomes of the bcf files, using the --regions. The rule i have created below works and does create a filtered bcf file.
rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} -o {output}"

However i would like to add the --regions parameter to the bcftools command. I have tried the shell part of the rule to several things
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions <chr1,chr2> -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {chr1,chr2} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1,chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1} {chr2}}
 {{chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions ['chr1','chr2'] -o {output}

But it doesn't work, and i am having trouble finding any suggestions to this specific issue.
I know it works if i just run bcftools filter in my command line, so using --regions should work.
Could anyone have any suggestions how to add --regions to the snakemake rule?

Have you tried expand() function? https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html?highlight=expand#the-expand-function

No i havent tried that, i didn't realize you could potentially use the expand() function for the shell part. But i will look into it

Something like the following would work:
rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    params:
        regions="chr1,chr2"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 --regions {params.regions} -v snps {input.bcf} -o {output}"

You could also include the regions as part of the output file name, rather than specifying it as a param.

Would a simple --regions \"chr1,chr2\" in the shell command not work?

You don't need to even quote it. However having it as a param allows you to set it elsewhere (e.g., the command line).


text:
Reading fasta sequence regions with faidx
Im coding in c++ and reading in different reference genomes to examine regions across the chromosomes a few hundred basepairs at a time.
To do it in c++ i use the library htslib and the command faidx_fetch_seq() which is similar to read specific regions using samtools faidx.
The problem is when im reading in the reference genomes, there are of course multiple regions consisting only of NNNNN (and so on). Since i only need the actual DNA sequence i wish to read in the entire file but removing the regions with NNNN.
I have tried reading in the regions in c++ and then just creating an if-statement to check the actual bases consist of N. This works but is incredibly slow.
So my question is: do anyone know how to only read in the actual DNA sequences, removing the NNN regions? either with samtools or with the c++ funciton faidx_fetch_seq().

Can you put "incredibly slow" in context? Slow is different for everyone, are we talking seconds, minutes, hours? Can you post some of your code that you wrote for this?

I just tested with a tiny program on a 400Mbp genome, reading in all the sequence while at the same time removing Ns. Took ~10.8 seconds, is that too slow? There are certainly ways to optimize this but you will be mostly bound by disk speed at this point, checking for Ns in a string of characters is not your bottleneck

Perhaps they are working with plant genomes, some of which will be unlikely to fit into the system memory of all but the most decked-out systems.

If you can't read in the entire sequence for a given reference genome:
(1) Use a script to generate a BED file of masked regions of low complexity (i.e. where you find Ns in your FASTA file). (2) Use bedops --complement with this mask file and the output of UCSC fetchChromSizes to generate a BED file containing regions that are unmasked. (3) Repeat for other reference genomes, as needed.
You can now loop over the unmasked regions in this file, making calls to samtools faidx or faidx_fetch_seq(). These calls will return unmasked sequence.
For performance or convenience, you could read chunks of sequence into a std::deque or std::vector buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.


text:
How to print alternative allele from VCF with htslib c++
I would like to print the alternative allele from a vcf file which i load in with htslib.
The description of some functions:
https://github.com/samtools/htslib/blob/develop/htslib/vcf.h
I have used bcf1_t objects by creating a header and reading the information.
htsFile *test_vcf = NULL;

  // creates header
  bcf_hdr_t *test_header = NULL;

  // initialize and allocate bcf1_t object
  bcf1_t *test_record = bcf_init();

  test_vcf = vcf_open("Testfile.vcf", "r");

  // returning a bcf_hdr_t struct
  test_header = bcf_hdr_read(test_vcf);

  while(bcf_read(test_vcf, test_header, test_record) == 0){
    std::cout << "chr" << bcf_hdr_id2name(test_header, test_record->rid) << ":" <<test_record->pos+1 << std::endl;
  }
  bcf_hdr_destroy(test_header);
  bcf_destroy(test_record);
  bcf_close(test_vcf);

But im not able to extract the reference or alternative allele information using this structure bcf1_t.
As far as i can read its perhaps possible with the structure and function:
bcf_dec_t
bcf_unpack()

But i am not able to get it to work since all the other functions seems limited to the bcf1_t struct.
Can anyone help with extracting the REF and ALT sequence?
Thanks for your time

see the code in https://github.com/samtools/htslib/blob/develop/vcf.c#L3104
 bcf_unpack((bcf1_t*)v, BCF_UN_ALL);// if you need genotypes too, might use another constant to avoid decoding all the structure
 (...)
 kputc('\t', s); // ALT
 if (v->n_allele > 1) {
        for (i = 1; i < v->n_allele; ++i) {
            if (i > 1) kputc(',', s);
            kputs(v->d.allele[i], s);
        }
    }

Thanks for reply, i'll try that!.

Okay, so as far as i understand based on this function, s is a kstring_t* which you can create using mem "kstring_t s = (kstring_t)&test_header->mem;", and then all the information is added to the pointer "s". But how are we then able to print or access the information in "s" ? I can see its a pointer declared with a reference, but then i should be able to print the values with *s. Again thanks for your time.

kstring_t is just the implementation of a dynamic string. you could just use cout << v->d.allele[i] << endl;

okay, that was my initial try, i just looked into it, and the problem is my n_allele for some reason is 0. So ill play around with it some more. Thanks again.


text:
General question: How to assign my SNP (Wheat, Illumina 15K) to specific chromosomes?
I have 2000 lines of wheat genotyped with a 15K illumina chip.
I also have the name of the SNPs.
I would like to know how could I find out which SNPs are in each chromosomes.

Answer from @devon-ryan, converted from comments:
Look in Biomart to see if there's SNP rsids for wheat.


text:
How to run Cox proportional Hazard model for each group in R?
My population is divided into 10 groups. I want to run the cox-proportional hazard model separately for each of these 10 groups so that in the end i can compare which groups has highest hazard ratio depending uopn the covariates.
How can i do this in R? Earlier i was using this formula.
Model_1 <- coxph(serve_object ~ GROUP + age+ Smoking+ Diabetes, data = DF
summary(Model_1)

But the output is very confusing. it gives hazard ratio of 9 out of 10 groups. I think it dealt group as covariate.

The fact that you supplied GROUP as a covariate, you are getting 9 coefficients as one of the ten levels is used as the reference level.
If you would like to fit one model per group (and in this scenario you would not be correcting for GROUP), you should first split your data frame and fit your models on the resulting list's elements:
lapply(split(DF, DF$GROUP),
       function(x) coxph(serve_object ~ age + Smoking + Diabetes, data = DF))

Example with the mtcars data set when grouped by number of cylinders, cyl having values of 4, 6 or 8:
> lapply(split(mtcars, mtcars$cyl),
+        function(x) lm(disp ~ hp, data = x))
$`4`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
    59.0369       0.5579

$`6`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   290.9178      -0.8799

$`8`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   320.2083       0.1572

Can you help me in interpreting result.
.$Group3
Call:
coxph(formula = so_ukb ~ age_at_recruitment.210220.0 + Smoking +
Hypertention, data = UKB_FILTRD)
                            coef exp(coef) se(coef)      z      p

age_at_recruitment.210220.0  0.06961   1.07208  0.00141  49.35 <2e-16
Smoking.L                   -0.21220   0.80880  0.01335 -15.90 <2e-16
Hypertention.L              -0.44221   0.64262  0.01361 -32.49 <2e-16
Likelihood ratio test=5142  on 3 df, p=< 2.2e-16
n= 448287, number of events= 11377
(3270 observations deleted due to missingness)

I am just confused that how to get one cumulative hazard ratio by adjusting model for these covariates. It is giving HR for age, smoking, hypertention separatively.. I want to see just cumulative.

Thank you for above lapply code, it worked so well.

Let's take the coefficient you get for age: The exp() of this coefficient is the change in in HR when age changes by one unit and while all the other variables are accounted for. In that sense, it takes all the variables into account so can be said to be "cumulative", although I am not sure such terminology exists for linear models. You might want to use interaction terms if your covariates depend on each other.

Can you give an idea of how we can export the summary(coxph) into csv?

or may be first convert it into data frame and then export df intom csv.

You can give a try with the broom package. Its tidy() function would convert output of lm() and many other related functions to data frames. Would probably work with the coxph()'s output.


text:
After completing the evaluation, how to predict PTM lysine or any other modification with unlabeled data for making a web server?
I have created a model with SVM classifier and evaluated it with 5-fold cross-validation. The dataset was sequence data containing lysine modified sites. Now I want to build a web server where anyone can input any sequence in FASTA format. Apparently, I have no label information at this point. So how can I will be able to extract features or use any feature selection method which requires label information? How can I extract or create the test set to predict? Suppose I have used f-score feature selection method in my model and selected 50 features from the given dataset. But for web server, having no label information, how can I select 50 features from the sequence that are provided by an user and then test it with my previously trained model?
As being a novice in the world of bioinformatics, I am seeking the answer which might be improper for others. But it will be very much helpful for me.
I am giving some link of some papers and its corresponding web servers.

iPTM-mLys: identifying multiple lysine PTM sites and their different type (https://academic.oup.com/bioinformatics/article/32/20/3116/2196525)
Web server link : http://www.jci-bioinfo.cn/iPTM-mLys
iHyd-PseCp: Identify hydroxyproline and hydroxylysine in proteins by incorporating sequence-coupled effects into general PseAAC
(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5190098/)
Web server link : http://www.jci-bioinfo.cn/iHyd-PseCp

There are too many questions really and not enough background, but I presume the gist of the question is this, you've a trained model and you want to know how to turn it into a website. So really you are looking for basic information on how to reproduce your model via a website. Is that correct?

Yes, it is @Michael. Basically I am using sequence-coupling effect to extract features and ANOVA F-test to select features where label information is necessary. Therefore, to reproduce my model in a website, I have to extract the query protein sequence given by an user and apply F-test again on it and later use it for detecting the class. Is this the process or there is other rules for it?
I have given two links that contain some background I hope.

Interesting application. I would simply suggest using Flask, which is described in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Obviously this is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution (same chapter). We're in lockdown so I don't have this book with me.

Thanks @Michael for suggesting me a path! Hope this will work. I look forward to reading it.

Interesting application. I would suggest using Python's Flask system, which is described in a few pages in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Django is not really recommended (don't say that on Stackoverflow however). O'Reilly do an entire book on Flask, but thats to make the site pretty.
Obviously, the explanation here is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution which is easy (same chapter). We're in lockdown so I don't have this book with me.
The difficult part will be web architecture because you shouldn't have the CPU/GPU on the same server as the quiry interface, otherwise you are prone to a denial of service attack or just sytem overload. To do this you need to consider  Representational state transfer (RESTful) and that is not easy at all (at least for me). Flask facilitates RESTful, but thats like systems admin stuff. I guess before getting into RESTful I'd see how popular your website is, but if its the same server as your institute website you'd need to be really careful because an external user could bring your institutional website down (you wouldn't be popular).

That would be very helpful. As I have to go through the book you have mentioned, I am not quite sure it is related to feature extraction or feature selection procedures. I just wanted to know that without any class labels, how can I classify a protein sequence. Such as an F-test requires class labels for statistical analysis. A query protein which is completely new and not included in my trained model's dataset, is to be classified as either positive or negative.


text:
Why my Seurat object has gene name as "gene-number"
I tried to analyze single-cell RNA sequencing data from PMID: 32849601. I downloaded fastq files. I also made a canine genome reference from https://www.ncbi.nlm.nih.gov/assembly/GCF_000002285.3/ according to a guideline from cellranger (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references#mkref-intro). I ran cellranger count as usual and created Seurat objects using features, barcodes, and matrix files based on the Seurat tutorial without modification.
I found all gene name has the suffix gene "-4" as below. [![enter image description here][1]][1].
Which step should I correct to have a plain gene name without the suffix? Since I have never found this issue before, it's really interesting. Is this related to the step of making a reference genome?
[1]: https://i.sstatic.net/oF1oi.png


text:
How to specifically select cluster on Featurescatter on Seurat?
Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger.
FYI, The filtered_feature_bc_matrix includes an aggregated population of tumor-infiltrating immune cells that are derived from 6 different syngeneic mouse tumors (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861). I downloaded fastq files and ran "cellranger count" and "cellranger aggr" commands to obtain combined data.
My strategy includes the use of the "feature scatter" function deposited in the library(Seurat). Finally, I ended up getting results below.

My question is that can we select a specific cluster(s) or identity on this plot? For example, I want to see correlation values in cluster 2 only.
The usage of Feature scatter is as follows.
FeatureScatter(
  object,
  feature1,
  feature2,
  cells = NULL,
  group.by = NULL,
  cols = NULL,
  pt.size = 1,
  shape.by = NULL,
  span = NULL,
  smooth = FALSE,
  combine = TRUE,
  slot = "data"
)
Do you know how to get the Pearson correlation of my two genes using Seurat? This is the best I can do.
Thank you.

You can use the cells argument to only display the data for a specific subset of cells.
For example to extract the values for a cluster 2:
# The cluster identity is stored in the meta.data
# In this example there is a clustering column called "RNA_snn_res.0.7"
# but any metadata column can be used

clust2.cells <- rownames(Seurat.data@meta.data[Seurat.data$RNA_snn_res.0.7 == 2, ])

FeatureScatter(Seurat.data, feature1 = "gene_A", feature2 = "gene_B", cells = clust2.cells)

Thank you so much.

You can also subset your Seurat object within the call to FeatureScatter, such as FeatureScatter(subset(Seurat.data, idents = "2"), "gene_A", "gene_B")


text:
How to concatename file names in Snakemake?
I want to pass a series of files (e.g., a, b) to a program prog like: prog a,b.
Is there a way to do this in Snakemake?
My attempt is:
rule run_prog:
  input: "a", "b"
  output: "output"
  run: "prog {input}"

But the actual command Snakemake generated is prog a b instead of prog a,b.
Any solution to pass a comma concatenated file names to Snakemake?

run: "prog {input[0]},{input[1]}" ?

@Chris_Rands that works, is there a general solution, as there can be ~20 input files. Thanks!

You can use join to convert the list of input file into a string.
I think this would work as expected.:
rule run_prog:
    input: "a", "b"
    output: "output"
    run: f"prog {','.join(input)}"

This works great. It works for f"prog {','.join(input.variable)}"  if I define variable under input: section. Can you explain why have to use input[:] instead of input?

Right, using input by itself is the same ! I edited my answer.


text:
Phasing problem using "samtools phase"
so I'm having a problem using samtools phase. I'm trying to phase my .bam archives that where generated by aligning my samples with my reference genome. Following the documentation I tried the command:
samtools phase -A  my_referenceGenome.fasta   my_sample_toBePhased.bam
The -A flag stands for dropping reads with ambiguous phase. The error I get is: [W::sam_read1] Parse error at line 1. I generated .bam files with and without the SAM header but the error persists.
I feel like I am missing something basic but that isn't clearly stated anywere. Does the my_referenceGenome.fasta archive needs to be previously phased? If so, how can I do it?
Thank you so much!

This is what I see in my samtools phase (v1.10)
Usage:   samtools phase [options] <in.bam>

samtools phase only takes a bam file as an input. But you are specifying a .fasta file as an input. If you want to specify the reference genome, you should be using the argument "--reference" as described here
      --reference FILE
           Reference sequence FASTA FILE [null]

Oh, you're right! I also tried that but it didn't work. I managed to make it work and I'll add an answer shortly. But before that I want to ask a little question:  in my pipeline I also remove pcr duplicates with the fixmate/markdup steps. Is it better to phase the reads before or after removing the duplicates?

I would recommend removing duplicates first.


text:
Running snakemake with multiple conda environments without rebuilding environments every time
I have some snakemake pipelines that have rules which run in different environments. Is it possible to run the pipeline without having to rebuild the environments from scratch every time when run in a different directory?
This is the way I would normally run snakemake:
snakemake --directory <my_output_dir> --use-conda
My snakefile might have rules that need to use different environments:
rule A:
   ....
   conda: envA.yml
   ....

rule B:
   ....
   conda: envB.yml
   ....

So if I wanted to do two different analyses with the same pipeline, is there anyway that I can point them to the same environments without having to install everything separately in each output directory?
The snakemake documentation makes it seem like they don't currently support this. You can either use 1) --conda-prefix when running snakemake and then point to a single env or 2) use the --use-conda directive and the env.yml files by rule, in which case it will create new envs any time you run it in a new output directory. So my case doesn't seem to be covered, but does anyone have a workaround to suggest? I often find conflicts between some of the R packages I use and python packages, so I tend to keep them in their own environments to avoid headache.
I'm starting to get neck-deep in conda environments at this point which is taking up space from my disk quota, but also I am bothered by the fact that I may actually want to run analysis of data from different projects with the same environments. If the environments are rebuilt at each execution, I do not know that they are the same, short of specifying the exact version info for everything I use in the env.yaml files (and even then, they could be slightly different if the repo has changed).
I have looked at this similar question: Running Snakemake in one single conda env
in which it is suggested there is a way to do such a thing, but I'm afraid I don't understand how to implement it without an example.
Thanks

Each YAML file will always create its own environment. The packages in each aren't downloaded and installed each time, if they're already installed elsewhere they're just hard-linked in. You'll want to use both --use-conda and --conda-prefix. The latter will point to the base directory in which all conda environments will be stored. This has the benefit that running the same workflow in different places won't result in a new env being created (by default snakemake will do so inside a .snakemake directory otherwise).

Wow, thank you. I had no idea that these two flags could be used together (I thought --conda-prefix could only be used with a named pre-built single environment) -- that solves the problem. I see what you mean about the hard-linking but it still ends up looking quite messy if you have all these .snakemake directories everywhere and it also seems to slow down some conda commands (like conda env list), at least for me.

@vantom, while I am not very familiar with snakemake on a whole, I do have extensive knowledge writing workflows in WDL as well as CWL. The approach both of these languages take to overcome a similar issue that you are describing, is to use containerization through ie docker. This approach allows you to bundle the specific tools you need for a specific task, making different images for different tasks. There are several large repositories of images online ie: Docker Hub, or more biologically focused: Biocontainers
Looking through the documentation, it does appear that snakemake supports  running your rule from within the context of a docker image. Without knowing how your compute environment is setup, I am not sure if there is any additional configuration you would need, however the documentation seems to suggest that it should be as easy as describing the name of the container. As long as you can run docker, (or an alternative like singularity) you should be able to avoid the current headache you are facing.
rule NAME:
    input:
        "table.txt"
    output:
        "plots/myplot.pdf"
    containers:
        "docker://joseespinosa/docker-r-ggplot2"
    script:
        "scripts/plot-stuff.R"

Thank you for the tip! I'm not currently using Docker, but it's useful to know this for the future.

For reproducibility, you should always used a tagged version of a container, that is:
docker://joseespinosa/docker-r-ggplot2:v0.1
If you don't it will pull/build the :latest tag and you run the risk your pipeline could break on an update.
See the tags menu in DockerHub; e.g. https://hub.docker.com/r/joseespinosa/docker-r-ggplot2/tags


text:
Is there a way to get the code from "github.mit.edu"?
everyone.
I am a postdoc who just begins to analyze single-cell RNA-seq data.
Nowadays, I found a really interesting paper (https://doi.org/10.1016/j.celrep.2018.10.047, PMID: 30404002). So, I decided to study bioinformatics by drawing the same graphs presented in the paper during the corona outbreak season. However, it was turned out that I can't get any access to the website in which code is deposited. Here is the address for the website.

Do you guys know how to get the code to learn bioinformatics and make graphs?
Thank you.

It is very strange that they put their code somewhere cannot be easily accessed. If they just put it under a normal Github repo, there would not be such a problem. One way to get the code is simply to email the corresponding author. People are usually responsive if you are interested in their work.

Thank you. I just sent the e-mail to PI. I hope to get the code information.

It is not uncommon for Institutions to set-up their own private github / gitlab instances to have more control over their code and what they do with it. Furthermore, it is quite likely that this was a private repo and someone forgot to change it to public. It happens. As the saying goes, never attribute to malice what can be explained by incompetence (or absent-mindedness) :)


text:
Error: Insufficient values in manual scale. 24 needed but only 1 provided?
folks.
Can you give me idea about an issue that I encountered?
I ran single-cell RNA sequencing data using Suerat in R, and when I tried to draw violin plot, there is an error as below.

VlnPlot(tumors, "Cd79a", "Foxp3", "CD14")

Error: Insufficient values in manual scale. 24 needed but only 1 provided.
Run rlang::last_error() to see where the error occurred.
However, I can draw the violin plot for Cd79a. The problem happened to only several genes. I haven't solved this problem yet. Do you know what should we do?
Thank you.