text:
At what stage of a transcriptome assembly is it better to perform read contaminant filter?
I'm trying to assemble a bivalve transcriptome. Since bivalves are filter feeders, their transcriptomes tend to be highly contaminated by bacteria, algae and whatnot. Since I pooled several transcriptomes, I have a high amount of reads (>2B reads). I also have a not-so-high quality reference genome.
Now, I'm not sure when to perform a BLAST-based filter on these reads:

Do I filter reads with BLAST before assembling them?
Do I assemble the transcriptome and then filter out contaminant contigs?
Since I have a reference genome, do I simply not care about contaminant reads because they will not align to it?

Having this amount of reads, I'd be inclined to first filter out those reads that don't have a match against Metazoa, and then to perform a second filter after assembly in order to retain only those contigs that match against Lophotrochozoa.
If you know resources on how to manage read contamination, feel free to share them.

First, to answer your question about mapping to a low-quality reference:
1. Mapping
For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome assembly does not contain any scaffolds from contaminating species like bacteria, algae, et cetera.
2. Reasons to not BLAST your raw reads
There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in #4, then reassembling, may work better.

3. Using multiple individuals of highly heterozygous species in assembly.
Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.
4. A strategy to filter transcripts based on mapping and BLAST
Given #1, #2, and #3 above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.
I would try this:

Assemble the transcriptome of each individual separately.

This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.

Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.

Use somthing like cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta.

Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.

Use samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam to only keep the transcripts that did not mapped to the reference. -f 4 is the magic flag here that keeps things that didn't map.

Make a list of all the transcripts that did not map to the reference...

samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt

...and get only the transcripts from your fasta file that are potential contaminants.

seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta

BLAST all of the unmapped.fasta and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this known_contaminants.fasta
Map all reads from all individuals to the known_contaminants.fasta.
Get all the reads that map to known contaminating transcripts, known_contaminants.fasta using:

samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt

Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
Reassemble the filtered fastq files.

This strategy has the following benefits:

Does not require weeks/months of compute time to BLAST 2 billion reads.
Ensures that all reads from the reference genome are left untouched.
Ensures that all reads from known contaminants are removed.
Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
Leverages contaminants from different individuals to filter reads from all individuals.

Thanks for your answer, there are some interesting ideas. I want to use BLAST mainly to filter out those reads/transcript that do have a bacterial/algal hit - on my first assembly attempt, I ended up with >20kb long transcripts that were just bacterial genomes. Do allelic differences distort the graphs even with a genome-based assembly? If I assemble single individuals first, I'll have 35 individuals - not counting the pooled larvae transcriptomes, which in themselves contain 50k individuals.

Allelic differences also distort the graph with genome-based assembly, causing many misassemblies and fragmenting the output. May I ask what sort of analysis that you're doing downstream with the transcriptomes? It will only take a small amount of time to make a shell script to assemble 35 transcriptomes compared to the amount of time it will take to discover all the incomplete and incorrect transcripts if all the reads are lumped into one assembly, IMHO.
I have updated the answer to give more specific instructions and the reasoning behind each step. Hope it helps!

Wow, that's a very thorough answer! Thank you very much, it seems a very clever strategy! I want to use these transcriptomes for a variety of analyses, GSEA, 3'-UTR motifs and other.

Hi @Dunois - I think I was being a bit dramatic above, but yes, point #3 is specifically relevant for de novo transcriptome assembly. You'll get more accurate results if a transcriptome is only assembled from reads from 1 individual. If you need different life stages or tissue types, my suggestion is to assemble the transcriptomes from different individuals independently, then remove duplicate transcripts using orthofinder. With that being said, if you have multiple libraries from different tissues from the same individual, then assemble all those together.


text:
Why DoHeatmap Does not show all genes in genes.use?
I am heatmaping a list of genes by DoHeatmap function in Seurat R package. I am sure I have 212 genes but heat map shows only a few of my genes
> DoHeatmap(
+        object = seurat,
+        genes.use = genes,
+        slim.col.label = TRUE,
+        remove.key = TRUE)
> dim(as.matrix(seurat@data))
[1] 12293   209

> length(genes)
[1] 212
> class(genes)
[1] "character"

> class(seurat)
[1] "seurat"
attr(,"package")
[1] "Seurat"
>

Please look at heat map, I was trying to plot 212 genes but you seeing only few genes are among my 212 genes
Please somebody save me from this confusion

This is probably a very naive question, but do more genes appear if you make the image larger? I mean, if you extend the image before exporting it, while it is being displayed in R?

Actually no, because before exporting from R looks the same :(

Yes, but what if you resize it before exporting? Does making the image bigger make more genes appear? I just want to rule out the obvious here.

Thank you, I just change the size of image before exposing to 2500*1991 but only rows became smaller and the order and number of genes are the same

Oh well. It was not very likely to be that simple.

https://github.com/satijalab/seurat/issues/287
In this post there is a answer for my question, however when I set scale to FALSE again the number of genes did not change but the shape of plot dramatically changed. So, maybe the issue is because of scaling that did not solve my problem

source: Stack Exchange

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

sort: votes


text:
Why is cell count stated relative to animal size?
Why do some single-cell RNA-seq papers give the number of cells sampled as a fraction of the size of the organ, animal, or embryo, as opposed to just saying how many cells they sampled?
For instance, Wagner et al.'s paper "Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo" says "For different developmental stages, we sampled 0.17x to 0.97x of the total cells per embryo, sufficient to detect cell states as rare as 0.1-0.5% of all cells." A pollster would never write "We sampled 4 out of every 10,000 registered voters in the USA"; they'd just give the number. Statistical results on cluster detection (e.g. Disentangling Gaussians by Adam Tauman Kalai, Ankur Moitra, and Gregory Valiant) are also phrased in terms of number of samples, not this ratio of sample size to organism size.
It's worth mentioning that these studies use low-yield techniques. In order to sequence 0.97 times the number of cells in a zebrafish embryo, Wagner et al. had to sequence at least 50-100 embryos. So, it's more like sampling 97 of 100 cells with replacement and less like sampling 97 of 100 cells without replacement.

The typical goal of scRNA-seq is to assess the possible cell-types and states that exist in an organism/tissue/etc. Given that, it can be useful to know what portion of the possible cells have been assayed, since that sets bounds on the number and size of cell populations not assayed. If the goal of a study is to assess the variability in a single cell-type, then you're correct that mentioning "0.8% of the population was sequenced" isn't actually informative. Note that "in a single cell-type" is important there, since if one is assaying multiple cell-types then it's important that enough of each are assayed given their prevalence in the relevant population.
Practically speaking, numbers like this are only ever written in papers if they look impressive, so no one is going to write about sequencing 0.000001% of the cells in an organism, since the reader/reviewer will then start to question the relevance of the results.

Could I ask you for more specifics on one point? You say this proportion they give, e.g. the 0.97, "sets bounds on the number and size of cell populations not assayed." How would you compute those bounds?

Normally you wouldn't compute those bounds, but if you really wanted to you'd have to assume some minimum size of a cell population and divide the total number of cells by it. Any numbers you'd get would be completely hand-wavy, of course, but that's how it goes in biology.

When you say "the total number of cells", you mean the total number of cells sequenced, right? If so, then this doesn't answer my question. My question is about why it is useful to represent the amount of cells as a fraction of the size of the system under study.

No, not the total number sequenced, the total number in the organism/tissue/etc.

It makes sense that you could use the total number in the organism/tissue to do what you're saying, but that is also not the ratio that I am asking about. I still don't see why Wagner et al.'s 0.97 is useful.

It's one of the ratios you're asking about. The 0.97 value makes sense if the fish were essentially genetically identical and from the same brood and the cells sampled from a given embryo are truly randomly selected.

I'm sorry, I can't understand your answer at all and I am not sure how to ask more clearly. Thanks for trying.

Your difficulty is in thinking that the discussion you linked to on cross-validated is relevant to the typical goals of scRNA-seq. It's not. Tissues/organisms aren't a population that we want to estimate some mean and confidence interval from, that's implicitly what scRNA-seq is trying not to do.

That's very far from explaining why the 0.97 is useful, but it's a good point and I may edit my question to make reference to some more relevant statistical discussion.

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


text:
How many reads do I need for hybrid assembly
I have Illumina and PacBio reads and I would like to use dbg2olc for hybrid assembly.
Part of dbg2olc is SelectLongestReads which select reads that sum to total_length of bases.
How many bases do I need for a 1.8G allotetraploid to remain?

I have never used dbg2olc. However other OLC assemblers managed to do a solid assembly from 20x of corrected long reads. But obviously, the more the better.
In you case of allotetraploid, you need to figure out if the heterozygosity is sufficient to separate haplotypes of not and adjust parameters of the assembler. For instance in bdelloid rotifers the tetraploid genome has an AA'BB' structure, all the assemblies were able to separate A and B, but not A and A' and B and B'. I guess an allotetraploid could have a similar structure.
It would be good to know what you are dealing with and what is your goal in the assembly. Do you want a single haplotype?
-- edit --
You can check the genome structure using short reads by smudgeplot, a tool I made (yes, it's a shameless plug).


text:
How I can test my hypothesises computationally
I have single cell RNA-seq data on about 2000 cells in 9 time point. I have clustered my cells in each time point by Seurat. I am seeing in some time points I have 3 clusters while in another time points I have 2 clusters of cells (please look at the plot).
One hypothesis would be that, a new cell type is being emerged at last time point. But maybe if we had another further time point this third cluster being disappeared because that has been arose just because of undifferentiated cells from the later time point (second hypothesis).
Another hypothesis would be this new third cluster is a new cells type that only temporarily expresses a set of genes but if we had cells for another further time point, we would not see the expression of these genes anymore and we had again 2 clusters. I am really struggling with this how I could test my hypothesises.
The final hypothesis would be that the third cluster is really a new cell type born from old cell types during the differentiation.
Like UAD algorithm addresses my question. I just tried FateID algorithm on my data but got stopped by
> pr <- prcurve(y,fb,dr,k=2,m="tsne",trthr=0.4,start=2)
Error in svd(xstar) : a dimension is zero

Any suggestion would be highly appreciated.

The answer is behind these papers:
The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution
Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo
Also this.

Please avoid link only answers, provide summary text from the links.

The papers are too difficult to follow but I just resealized that finding parent child clusters likely would work.

This is what packages like FateID and Monocle are for, namely taking single-cell RNAseq data and inferring differentiation trajectories from it. Don't try to reinvent the wheel on this, there are a number of packages out there to do this sort of thing and they're going to have pretty complicated methods that you're really not going to want to reinvent unless you absolutely have to.


text:
Get results of keyword search on Pfam via python script
I'm interested in all proteins that are in any way associated with Danio rerio. I decided to look them up at Pfam data base and when I just make a keyword search, I get a a nice list which looks like this http://pfam.xfam.org/search/keyword?query=danio+rerio&submit=Submit.
Important thing to note is that most of the associations of these protein families can be found in seq_info column.
My question is: how to extract that list from a python script? Moreover, I want to make a script that when run, will fetch the latest release, filter out the proteins that cannot be found in D. rerio and do some additional processing.
I've been searching for ways to directly download files from Pfam ftp server and mine for the proteins that are associated with D. rerio, however, I didn't find a single way to do it. I simply cannot see where do they store the information about species in which the protein can be found.
I would like to retrieve everything about protein families of interest, because the next step I will need to do is to find the DNA binding proteins and for finding a way to filter out the ones of no interest, I would liketo have as much data as possible.

If I understand your question correctly, then you can use the Pfam proteom search to get the data you are looking for.
Go to the

Pfam website -> Browse -> select the letter D from proteoms category
  -> D. rerio

There you can find a download link under the section Domain composition:

You can download a list of all regions for this proteome from our FTP
  site.

The download link seems to be bound to the species ID and if you want to automate your data retrieval you may be able to simply download the file from there in your python script.

Exactly what I needed. Although, the species ID is not that informative but I think I can manage that easier in the future. Thank you!

One thing to note though is that seq_info at Pfam is retrieved from UniProt database so in order to retrieve that data, one needs to make a mapping from list of proteoms and uniprot.gz file at the Pfam ftp server. The tricky part comes if the two databases are not synced. But that's not an issue that should be discussed here.


text:
A strange error while clustering of single cell RNA-seq using URD package
I am using URD R package on my single cell RNA-seq data hoping to be able to trace gene expression changes during the development. After dimension reduction step by PCA and diffusion map whatever I am trying to cluster some cells, but whatever I do I am getting this error
> object <- graphClustering(object, num.nn = c(5,8,10,15), method="Louvain", do.jaccard = TRUE)
Error in graphClustering(object, num.nn = c(5, 8, 10, 15), method = "Louvain",  :
  unused arguments (num.nn = c(5, 8, 10, 15), method = "Louvain", do.jaccard = TRUE)

> print(object.6s)
URD object: 12451 genes x 1556 cells.

object.6s.mnn <- graphClustering(object.6s)
Show Traceback Rerun with Debug Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘graphClustering’ for signature ‘"URD"’

I tried to change Louvain, jaccard, ... but I am still getting this error. How can I get rid of it?

The error message says that you should remove the num.nn, method and do.jaccard argument. Did you try it?

Thank you, yes I removed all except for my URD object but the error is permanent :(

Can you post the new command and the errors? Could you add the 'print(object)' result?

Please [edit] the question to add the code (in code formatting) . It seems like the error is different now. I'll research a bit

The problem is they define a new class but they don't define methods for said class. The package would need (IMHO) more test and add methods for said class. It would be still possible to work with functions if they defined a generic class, but that's beyond what we can answer/solve here, they need to modify the package.

Also, please consider asking to the right people next time. If you post an issue in the package it signals that there is something wrong with the package(you followed the tutorial as you did and you found a problem), while if you post it here I understand that there is some problem while using the package (you are trying to do something not covered by the help file or the tutorial and don't know how to do it). Either way if it is a "simple" informatics error probably Bioinformatics.SE is not the best place to post it.

Thank you I also posted in GitHub but URD package developers did to reply

Yes, I saw, it but don't expect people to answer the emails or questions in 5 minutes. The developers (and people on the forums, here on biostars or in support.bioconductor) have other things to do and are answering you on their own time, because they want to help you. To help you they need to spend some time figuring the error and correcting it for you, so it takes some time to give you an answer. I hope you get your answers and solution to these problem.

definitely you are alright, I am also thankful for your time

This solved the error
object <- URD::graphClustering(object)

By this function I am telling R to use the graphClustering() function from the URD namespace. I probably had loaded another package with a different graphClustering() function - or defined one myself - which is taking precedence over URD's function.

It's nice that you have shared a solution to your problem. But since it was not very clear what is the actual problem, this answer is not really explaining anything.

I added the possible reason


text:
Are these standard species abbreviations and how to look up others?
Reactome uses a three letter code for each species which are listed below.
Are the letter codes used by Reactome a standard?
If so is there a web page or download file that can validate the code and translate it to a common name?
 HSA  Homo sapiens
 SSC  Sus scrofa
 GGA  Gallus gallus
 DME  Drosophila melanogaster
 RNO  Rattus norvegicus
 DDI  Dictyostelium discoideum
 MMU  Mus musculus
 ATH  Arabidopsis thaliana
 TGU  Taeniopygia guttata
 DRE  Danio rerio
 OSA  Oryza sativa
 BTA  Bos taurus
 XTR  Xenopus tropicalis
 CFA  Canis familiaris
 SCE  Saccharomyces cerevisiae
 PFA  Plasmodium falciparum
 CEL  Caenorhabditis elegans
 SPO  Schizosaccharomyces pombe
 ROT  Rotavirus
 FLU  Influenza A virus
 HCV  Hepatitis C Virus
 HBV  Hepatitis B virus
 HER  Human herpesvirus 1
 NME  Neisseria meningitidis serogroup B
 NME  Neisseria meningitidis
 CTR  Chlamydia trachomatis
 STY  Salmonella typhimurium
 STY  Salmonella enterica
 ECO  Escherichia coli
 OCU  Oryctolagus cuniculus
 XLA  Xenopus laevis
 CAE  Cercopithecus aethiops
 CGR  Cricetulus griseus
 HER  Human herpesvirus 8
 MCV  Molluscum contagiosum virus subtype 1
 MCV  Molluscum contagiosum virus
 HCY  Human cytomegalovirus
 HIV  Human immunodeficiency virus 1
 PCH  Penicillium chrysogenum
 TAE  Triticum aestivum
 VRA  Vigna radiata var. radiata
 VRA  Vigna radiata
 CFS  Crithidia fasciculata
 FCA  Felis catus
 AMA  Arenicola marina
 MGA  Meleagris gallopavo
 NGO  Neisseria gonorrhoeae
 SAU  Staphylococcus aureus
 MVI  Measles virus
 HAL  Human alphaherpesvirus 2

EDIT
Not a complete answer as this only covers a subset of all organisms but at present better than nothing.
Using KEGG Organisms
https://www.genome.jp/kegg-bin/get_htext
Enter letter code, e.g. GGA

EDIT - After reading answer by conchoecia
The graph database schema for Reactome only has two labels that use the 3 letter code which is named abbreviation in the schema:
Species and ReferenceTherapeutic
Species has taxid  (Taxon Identifier)
abbreviation     taxId       displayName
AMA              6344        Arenicola marina
ATH              3702        Arabidopsis thaliana
BAN              1392        Bacillus anthracis
BTA              9913        Bos taurus
CAE              9534        Cercopithecus aethiops
CAL              5476        Candida albicans
CBO              1491        Clostridium botulinum
CEL              6239        Caenorhabditis elegans
CFA              9615        Canis familiaris
CFS              5656        Crithidia fasciculata
CGR              10029       Cricetulus griseus
CPH              10703       Corynephage beta
CPO              10141       Cavia porcellus
CSB              60711       Chlorocebus sabaeus
CTE              1513        Clostridium tetani
CTR              813         Chlamydia trachomatis
CVI              10243       Cowpox virus
DDI              44689       Dictyostelium discoideum
DME              7227        Drosophila melanogaster
DRE              7955        Danio rerio
ECO              562         Escherichia coli
FCA              9685        Felis catus
FLU              11320       Influenza A virus
GGA              9031        Gallus gallus
HAL              10310       Human alphaherpesvirus 2
HAM              6706        Homarus americanus
HBV              10407       Hepatitis B virus
HCV              11103       Hepatitis C Virus
HCV              31649       Hepatitis C virus genotype 2a
HCV              31646       Hepatitis C virus subtype 1a
HCY              10359       Human cytomegalovirus
HER              10298       Human herpesvirus 1
HER              37296       Human herpesvirus 8
HGA              10376       Human gammaherpesvirus 4
HIV              11676       Human immunodeficiency virus 1
HPV              333760      Human papillomavirus type 16
HPV              337041      Alphapapillomavirus 9
HSA              9606        Homo sapiens
IBR              11120       Infectious bronchitis virus
LMO              1639        Listeria monocytogenes
LMO              1906951     Listeria monocytogenes serovar 1/2a
LPN              446         Legionella pneumophila
MCV              10280       Molluscum contagiosum virus subtype 1
MCV              10279       Molluscum contagiosum virus
MGA              9013        Meleagris gallopavo
MML              9544        Macaca mulatta
MMU              10090       Mus musculus
MTU              1773        Mycobacterium tuberculosis
MTU              83332       Mycobacterium tuberculosis H37Rv
MVI              11234       Measles virus
NGO              485         Neisseria gonorrhoeae
NME              491         Neisseria meningitidis serogroup B
NME              487         Neisseria meningitidis
OAR              9940        Ovis aries
OCU              9986        Oryctolagus cuniculus
OSA              4530        Oryza sativa
PCH              5076        Penicillium chrysogenum
PFA              5833        Plasmodium falciparum
RNO              10116       Rattus norvegicus
ROT              10912       Rotavirus
SAU              1280        Staphylococcus aureus
SCE              4932        Saccharomyces cerevisiae
SEV              11191       Sendai virus
SPO              4896        Schizosaccharomyces pombe
SSC              9823        Sus scrofa
STY              90371       Salmonella typhimurium
STY              28901       Salmonella enterica
TAE              4565        Triticum aestivum
TGO              5811        Toxoplasma gondii
TGU              59729       Taeniopygia guttata
VRA              3916        Vigna radiata var. radiata
VRA              157791      Vigna radiata
VVI              10245       Vaccinia virus
XLA              8355        Xenopus laevis
XTR              8364        Xenopus tropicalis

While ReferenceTherapeutic does not.
The species can be identified using the NCBI taxon identifier, e.g.
https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=6344

EDIT - After reading answer by user172818
Searching the file https://www.uniprot.org/docs/speclist.txt for 6344
AREMA E    6344: N=Arenicola marina
                 C=Lugworm
                 S=Lumbricus marinus

or using the URL: https://www.uniprot.org/taxonomy/6344

and with reference to

I think databases should stop inventing new species abbreviations.

The Wikipedia article on Lugworm shows the many Taxon identifiers in multiple databases.

No, there are no standardized three-letter abbreviations for species names. Three-letter abbreviations with roman letters is 26^3 combinations, which is only 17,576 possible combinations. As a result there are not enough three-letter abbreviations for all living species.
All three-letter abbreviations are either specific to a database or perhaps are shared between a database or two. There is currently no substitute for a binomial name (although NCBI taxid is a convenient parsable substitute for the species included in NCBI to date).

These are not. The closest to "standard" is the 5-character abbreviation by Swissprot. It names 25,886 species and has been used for decades. It is also easy to remember for some common species such as HUMAN, RAT, MOUSE and HORSE. I think databases should stop inventing new species abbreviations.


text:
Finding gene length using ensembl ID
I want to find the length of a list of genes, of Homo sapiens, that is reported in the GEO database. I have gathered the ensembl id's of those genes.
I understand this information can be parsed from the start and end position given in GTF file available in ensembl.org.
I would like to ask for advice on how to programmatically download the gtf files with ensembl id's.
On the ensembl webpage a couple of options are listed, API and ftp access. I would like to know if sample scripts are available.
Update: I found a script here, it shows how to parse the data from gtf file. Could someone help in understanding how the code can be modified to obtain gene length from start and end position?
Edit 1:
By start and end position, I refer to the start and end values given in the "fields" section here and the file from which I want to parse the values is the second last file available in this link.
For instance, the end - start value of the gene with ensembl id ENSG00000131482,
gene :G6PC
Edit 2: I get the gene names and the expression values from the SOFT files documented in GEO, link
The following is a sample output from SOFT file. It gives me the probe_set(ID_REF) , Gene name and expression values from samples (GSMXXXXXX). Here we see fro G6PC, values are reported from two probe_sets. Can we get the transcript id's from the probe_sets?
Obtaining the transcript id will help me in finding the length of the specific transcript.
> Table(gds)[2449,0:5]
           ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
2449 1555612_s_at       G6PC   4.39472   4.52235   4.26869
> Table(gds)[16399,0:5]
         ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
16399 206952_at       G6PC   2.76693   3.45107    2.9815

What do you mean by gene length? Is that the distance from the start of the first CDS to the end of the last CDS in the gene, the length of the longest transcript, or something else?

The start and the end position given in gtf files(I'm trying to obtain thi length), does it include both the coding domain and the non-coding domain? I'm a beginner in this field, excuse me for the naive questions

That would then be coding, non-coding and intronic regions. What do you need these values for?

Can you give us an example? Show us a few lines of the GTF file and the output you expect to see from them.

I'm trying to see whether there is a correlation between gene length and the gene expression values reported from micro array studies performed in different platforms.

OK, but please [edit] your question and show us an example of your input data so we can test our solutions. You also still haven't clarified if you mean the length of the gene (which might include regulatory regions), the length of each transcript (different transcripts will be differentially expressed) or the length of the CDS of each transcript.

@terdon Please check my edit. I don't have the ensembl id's of the transcript. From the GDS files available from the GEO webpage, I have parsed the gene list. I could find only the gene names from those files and from what I understand those files are not annotated with transcript id's.

There are 4 different gtf files there, each of which has different types of data. I still haven't understood if you really want the start/end position of the gene or of each transcript. If you're looking at expression, you probably need the transcript but you keep referring to genes. You might want to come into our chat room so we can discuss this and clarify your question. If you do, please ping me (write @terdon).

I can save you time and say that microarray signal is unlikely to be affected by gene length.

@DevonRyan Realized that after doing the whole task. Thanks anyway:)

If you really just want the start and end positions of all genes in a list of genes from Homo_sapiens.GRCh38.93.gtf.gz, you could do this in awk. For example, if your list of genes is:
$ cat genes
ENSG00000131482
ENSG00000223972
ENSG00000228794
ENSG00000187961

You can do:
$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$4,$5
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 11869 14409
ENSG00000228794 825138 859446
ENSG00000187961 960587 965715
ENSG00000131482 42900797 42913369

Explanation:

if(NR==FNR){ wantedGenes[$1]++ } : NR is the input line number and FNR is the current line number of this input file. The two are only equal while reading the first file. So, if this is the first file, save the first field of each line (the gene name) in the array wantedGenes.
if($3=="gene"){ : if this line is a 'gene' entry.
gsub(/";*/,"",$10); : the gene names in the GTF file are quoted and have an ; after the quote. For example: "ENSG00000187961";. Since we don't have this in the list, we need to remove it so we can match the gene name. This gsub command will remove all " followed by 0 or more ; (the ;* means 0 or more occurrences of ;) from the 10th field of each line. So, it will replace "ENSG00000187961"; with ENSG00000187961.
if($10 in wantedGenes){print $10,$4,$5} : if the 10th field (the modified gene name from above) is in the wantedGenes array, then print the 10th field (the gene name), the 4th (start) and the 5th field (end).

You can easily modify this to print the length instead of the start and end positions:
$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$5-$4
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 2540
ENSG00000228794 34308
ENSG00000187961 5128
ENSG00000131482 12572

As requested in chat, here's a python version of the same thing:
#!/usr/bin/env python3
import re
import sys

wantedGenes = []
with open(sys.argv[1], mode='r') as f:
    wantedGenes = f.readlines()

wantedGenes = [ x.rstrip() for x in wantedGenes]

print('# Gene\tstart\tend\tlength')
with open(sys.argv[2], mode='r') as f:
    for line in f:
        if not line.startswith('#'):
            fields = line.split('\t')
            annotations = fields[8].split(' ')
            geneName = re.sub('[";]', '', annotations[1])
            if geneName in wantedGenes and fields[2] == 'gene':
                print('%s\t%s\t%s\t%d' %
                      (geneName, fields[3],fields[4],
                       int( int(fields[4]) - int(fields[3]) )
                      ))

Save the file as foo.py and run with the file names as arguments:
python foo.py geneList.txt Homo_sapiens.GRCh38.93.gtf

You could also do this in R with the GenomicFeatures library.
library(GenomicFeatures)

## make TxDb from GTF file
txdb <- makeTxDbFromGFF('Homo_sapiens.GRCh38.93.gtf')

## get gene information
all.genes <- genes(txdb)

## import your list of gene names
my.genes <- c('ENSG00000141510','ENSG00000184571','ENSG00000011007')

## get the length of each of those genes
my.genes.lengths <- width(all.genes[my.genes])
## put the names back on the lengths
names(my.genes.lengths) <- my.genes

## print lengths
print(my.genes.lengths)

ENSG00000141510 ENSG00000184571 ENSG00000011007
      25772           55687           18905

You can do this directly on BioMart site from EBI as well. Using a list of ENSG ids as input, you can get the actual length of each individual transcript, as well as the length of the genomic region coding said transcripts. These are the attributes I used:
Gene stable ID
Transcript stable ID
Transcript length (including UTRs and CDS)
Gene start (bp)
Gene end (bp)

The output is a table with the relevant data:
Gene ID Transcript ID   Transcript length (including UTRs and CDS)  Gene start (bp) Gene end (bp)
ENSG00000146648 ENST00000455089 3844    55019021    55211628
ENSG00000146648 ENST00000342916 2239    55019021    55211628

You can change the build as well if you don't want to use GRCh38.

If you are an R user, then you can easily do so with biomaRt
library("biomaRt")
library("dplyr")

genes <- read.delim("genes.file.txt", header=T)

head(genes)
          gene_id
1 ENSG00000000003
2 ENSG00000000005
3 ENSG00000000419
4 ENSG00000000457
5 ENSG00000000460
6 ENSG00000000938

## Build a biomart query
# In the example below, I use the human gene annotation from Ensembl release 82 located on "sep2015.archive.ensembl.org", More about the ensembl_build can be found on "http://www.ensembl.org/info/website/archives/index.html"

dataset="hsapiens_gene_ensembl"
mart = biomaRt::useMart("ENSEMBL_MART_ENSEMBL", dataset = dataset,
host = paste0("sep2015", ".archive.ensembl.org"), path = "/biomart/martservice", archive = FALSE)

annotations <- biomaRt::getBM(mart = mart, attributes=c("ensembl_gene_id", "external_gene_name", "start_position", "end_position"))
annotations <- dplyr::transmute(annotations, ensembl_gene_id, external_gene_name, gene_length = end_position - start_position)

# Filter and re-order gene.annotations to match the order in your input genes list
final.genes <- annotations %>% dplyr::filter(ensembl_gene_id %in% genes$gene_id)
final.genes <- final.genes[order(match(final.genes$ensembl_gene_id, genes$gene_id)),]; rownames(final.genes) <-NULL

head(final.genes)
  ensembl_gene_id external_gene_name gene_length
1 ENSG00000000003             TSPAN6       12882
2 ENSG00000000005               TNMD       15083
3 ENSG00000000419               DPM1       23688
4 ENSG00000000457              SCYL3       44636
5 ENSG00000000460           C1orf112      192073
6 ENSG00000000938                FGR       23213


text:
Fitting a principal curve into a diffusion map
I returned a Seurat object with a diffusion map by
 seurat_object <- RunDiffusion(seurat_object,genes.use = seurat_object@var.genes)

as below

Now, how I could know the order of a cell's projection onto this principal curve representing its predicted progression through the development by fitting a principal curve on my diffusion map?
I know Seurat slot for DM1 and DM2 is inside the seurat@dr
by
fit=principal_curve(as.matrix(seurat_object@dr$dm@cell.embeddings))
plot(fit)
line(fit)
points(fit)

I draw this graph of my cells

But how I know the first and last cells on curve? I mean how I know cells from each cluster is being placed where on the curve? By fit$ord I will obtain orders of cells on the line but how I know the first and end of the line?????

If fit$ord gives you the order of cells on the curve, you should be able to obtain the first and last cells by doing the following:
first_cell_index <- which.min(fit$ord)
last_cell_index <- which.max(fit$ord)

However, I'm not sure this is so informative - if you're looking to get representative gene expression for the start and end of the curve, you should really be looking at the cluster of cells towards the start and end of the curve - for example, by comparing cells in clusters 1 and 3 that you showed in your first plot. A single cell will suffer from dropout and not tell you a lot about the cluster.

Thank you, fit$lambda gives me pseudotime amount; Now I want to have a violin plot for lambda for each cell and each cluster


text:
Counting letters in phylip alignment columns with Biopython
I have been using python 3.6 and biopython 1.72 to work with protein data files. I am using a protein sequence file (phylip format), for example:
    14    678
Zebrafish  LSSCGVVSGD LISVILPASS LEETQTSSAA AHQTHTDQQA GGSHVSSSSS
Fugu       LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS-
Chicken    LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS-
Zebra      LASYGIVPGD LICLLLEETD AKSSLPPPSS SSPPSLQNGR ETSTLTPKNQ

I want to manipulate the file in a way so that;

I may get the number of columns which have all the same letters (amino acids)
The number of columns which have at least 2 different amino acids (with no dashes).
The number of columns which have any single dash.

I have iterated over the columns of the file with:
alignment = AlignIO.read(open("example.phy"), "phylip")
def column(alignment):
    n = len(alignment[0].seq)

    for i in range(n):
        c = []

        for record in alignment:
            c.append(record.seq[i])

        yield c

This returns the columns of the file, but I don't understand how I can manipulate them further.
Can I get some helpful suggestions?

BTW alignment and sequence-alignment tags should probably be synonyms

I suggest you use zip() to transform your data into columns and then apply the following counts similar to heathobrien suggests:
from Bio import AlignIO

alignment = AlignIO.read(open('example.phy'), 'phylip')

var1 = var2 = var3 = 0
for col in zip(*alignment):
    num_unique = len(set(col))
    if num_unique == 1 and col[0] != '-':
        var1 += 1
    elif num_unique > 1 and '-' not in col:
        var2 += 1
    elif '-' in col:  # assumes 1 or more dashes
        var3 += 1

print('Number of columns with the same amino acid: {}\n'
      'Number of columns with at least 2 amino acids (no gaps): {}\n'
      'Number of columns with at least one gap: {}'
       .format(var1, var2, var3))

For your example data, this outputs:
Number of columns with the same amino acid: 9
Number of columns with at least 2 amino acids (no gaps): 40
Number of columns with at least one gap: 1

using something like for col in column(alignment):

I may get the number of columns which have all the same letters (amino acids)
len(set(col)) == 1

The number of columns which have at least 2 different amino acids (with no dashes):
len(set(col)) > 1 and '-' not in col

The number of columns which have any single dash:
'-' in col

or if you really want columns with a single dash:
    col.count('-') == 1

Earlier I have used a for loop to iterate over the columns for c in column(alignment): and get the columns with dashes if "-" in c:         list2.append(c)
and without dashes     else:         list3.append(c)
in two separate lists. But when it comes to count the total number of columns with dashes, it results out the number of dashes for each of the column in file and keeps on printing unless all the columns are iterated (because of for loop).

It would be helpful if you can elaborate me what you have suggested because I cannot follow up the idea. I still have my problem. 1) I may get the number of columns which have all the same letters (amino acids) **2)**The number of columns which have at least 2 different amino acids (with no dashes). **3)**The number of columns which have any single dash. –


text:
Getting all NNI trees of a parsimony tree
I have an aligned protein sequence file which I have been using for reconstructing a parsimonious tree. I am currently using NNITreeSearcher._get_neighbors method from Biopython 1.72 but it's way to find one best scored tree only.
def _nni(self, starting_tree, alignment):
    """Search for the best parsimony tree using the NNI algorithm (PRIVATE)."""
    best_tree = starting_tree
    while True:
        best_score = self.scorer.get_score(best_tree, alignment)
        temp = best_score
        for t in self._get_neighbors(best_tree):
            score = self.scorer.get_score(t, alignment)
            if score < best_score:
                best_score = score
                best_tree = t
        # stop if no smaller score exist
        if best_score >= temp:
            break
    return best_tree

I need to generate all NNI-neighbours of the trees (two trees are NNI-neighbours if one can be transformed into another by one nearest neighbour interchange operation).

I want to know how can I get total number of trees generated before the best tree is reached.
Is it possible to return all the trees generated before the best tree is reached, so that the all/required number of trees(other than the best tree) may be printed?

As I've studies the 'TreeConstruction' module of Biopython, I came to know that '_get_neighbors' method is the one that is used to get all neighbor trees of the given tree. This is all I've concluded from my findings. If anyone wants to add more, that'll be positively appreciated.


text:
Get a list of PDB ids from a list of sequence protein number?
I am working on a project about phosphorylation, and I am creating a database of some proteins. I have a list of protein accession ids for which I need to get the corresponding PDB file for each one to map it to PSM matrix. I want to get the PDB file using Python or BioPython tools. Any suggestions will be appreciated.

Thanks for your reply. I have a list of protein ids, and I want to retrieve PDB code for each one! I made the protein ids database with the phosphorylation cite for each position. However, I need PDB code because I want to map it to PSM matrix. I have no idea how to retrieve PDB code! Thanks again

Please include a few IDs and the expected outcome, i.e. ID: P26378 should map to PDB: 1FXL, etc. If you just need to do it a few times, you could use the UniProt web interface, just upload a bunch of identifiers and get the corresponding PDB IDs. https://www.uniprot.org/uploadlists/

Thanks man, I really appreciate your reply. I am currently using this website to map protein accession ids to  PDB. However, I want to do it by the code because I have a big list of protein ids. I want to get PDB by code so I can map it to PSSM to get the structure of the proteins. I want to do it in the one piece of code,not separately! Again Thanks for your kind and helpful reply!

A typical Uniprot URL would look like the example here:
https://www.uniprot.org/uniprot/?query=ID:P26378&format=tab&columns=id%2Cdatabase(PDB)
The response will be:

Entry Cross-reference (PDB)
P26378    1FXL;1G2E;

We want our query to include the protein ID, the output format should be tab separated and the only contain the columns id and database (PDB).
To get the PDB IDs you would need to split the response, take the last line, split this line after each tab and split it again after each semi-colon.
Using Python's Requests library the whole code would be.
import requests
import time

uniprot_ids = ['P26378', 'O35433', 'Q02910']
url = 'https://www.uniprot.org/uniprot/'

protein_to_pdb = {}
for protein in uniprot_ids:
    params = {
        'format': 'tab',
        'query': 'ID:{}'.format(protein),
        'columns': 'id,database(PDB)'
    }
    contact = ""  # Please set your email address here.
    headers = {'User-Agent': 'Python {}'.format(contact)}
    r = requests.get(url, params=params, headers=headers)

    protein_to_pdb[protein] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')
    protein_to_pdb[protein].pop(-1)
    time.sleep(1)  # be respectful and don't overwhelm the server with requests

print(protein_to_pdb)

To accelerate things, you could store your protein_to_pdb dict in a file and only retrieve the new ones
Checking if the returned ID matches your input ID is always a good idea

Thanks man, I really appreciate your answer and your time with me. Please I want to read a file which has a list of protein ids as df ,and retrieve all PDB codes for each one. I am trying to edit on your code by using :"df.iloc[i,1] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')" ,but there is an error which is "list index out of range". Please can you  help me with that?

I did it by using :"df.iloc[i,1] = r" instead of using a list,and I get a file with all protein ids and PDB codes. Is that correct? Or I am missing something! Thank you so much for your answer; it is really helping me a lot. I really appreciate it .

Try a more simple example, just copy a few ids in the list and see if it works. Without seeing your dataframe it is hard to guess where the problem occurs.

You could use mapping between UniProt and PDB entries from the SIFTS project.
The uniprot_pdb.csv file there has UniProt accession names followed by a semicolon-separated list of PDB four letter codes.


text:
Gathering data on bacterial organism growth conditions
let's say I have a list of organism names like the example below
Achromobacter xylosoxidans
Acidithiobacillus ferrooxidans
Chloroflexus aurantiacus
Clostridium perfringens
Aquaspirillum arcticum

The goal is to identify the optimal growth conditions for each organism--meaning the optimal growth temperature and pressure. Obviously I can just google this information if the list is composed of a few organisms or even few tens of organisms. However, assuming that the list is too large, let's say several 100s of names, what would be an efficient means to gathering this information.
In searching for this info, I have not come across a database that has the information that I'm looking for in one place. Instead I find that this information is scattered across multiple sites, databases, scientific publications, protocols, etc. For each organism, I basically have to search through and extract this manually, and I was wondering if anyone can suggest some way to automate this. I am familiar with Python3, ideally I would accomplish this with python.
I though about implementing some sort of a web scraping solution, where I provide the name of the organism, as well as a few key words and see what happens. However, many times this data is not well tabulated and the name and whenever you do find information on the organism, you need to read to be able to interpret the information to get your answer...
maybe I'm not looking in the right place, or my approach to this problem is wrong/backwards/etc, any suggestions will help.
I haven't been able to find a centralized place to gather this information, that's why I posted here. The closest thing I've found is the ATCC.

https://www.biostars.org/p/235851/

yuo need to parse all of their product sheet, colored pdfs .....https://www.atcc.org/api/pdf/product-sheet?id=49176

You could use the BacDive API, or you could download the kgx formatted files from KG Microbe.


text:
VlnPlot problem
I was confused about the VlnPlot that it shows like the screenshot. It only shows the dot but no distribution of the data, is that normal?

Hello Hua, welcome to bioinformatics.SE. If you search google images for "Vln Plot", you find how it should look like, for instance. We are happy to help with questions, but we kindly ask people to devote some effort in resolving problem by their own and showing the research effort in the question. Also note that we don't see your data, it's very hard to guess what is wrong with your plots without it. My guess is that you have vast vast majority of data points == 0.

But is it possible to show the distribution for lowly expressed genes? How can I manually set the parameters in vlnplot() function to show the contours?

I assume you are referring to VlnPlot() of Seurat. The reason that you are getting such a plot is because your distribution is highly skewed, most of your points are zeros, not surprising for scRNA-seq. And the plot is actually showing your distributions, it is just the points at zero stacking on top of each other.
You can double check this by plotting a histogram for each of the clusters.
An example violin plot:

Cluster 0, which is mostly zeros, kind of looks like yours.

Cluster 6, which I believe is what you would expect to see, contains a lot of cells that express the gene and much less skewed, hence the shape denoting the density.

So what you are seeing is normal, your gene is not expressed/not detected is all.

Try using ridgeplots instead of violin plots.  Ridge plots will show you that bump at zero


text:
Convert TPM-normalized matrices back to UMI in python
I want to process a plenty of scRNAseq datasets in python, and I want to run TMM normalization on all of them. Unfortunately some of them are already normalized with TPM method. Is there a way of either converting it back to UMImatrix, so that TMM normalization could be applied to it or could we just apply TMM directly to TPM-normalized matrices? Would that work? Any suggestions would be greatly appreciated.

It depends what your downstream requirement is. If you want to do count based statisitcs, then you definately shouldn't just do TMM normalisation on the TPM data as count statistics only work on real count data.
If you are intersted in doing something like clustering or dimensionality reduction, then I wouldn't use TMM, I'd probably use rLog, which includes normalisation.
You can have a crack at converting the TPMs back to something that looks like counts, but I not sure how well its going to work. To convert them back to REAL counts, you are going to need to know: the total number UMIs in the original dataset and the precise gene length used for each gene. Its posible the data set provides this ... Otherwise you are going to have problems.
$\frac{TPM_i}{1,000,000}$ is the fraction of all transcripts in a sample that belong to a given gene.
If $$fL_i = \frac{length_i}{\Sigma_{j\in G} length_j} $$
is the fraction of the transcribed genome in a gene i. Thus $fR_i = fL_i \times \frac{TPM_i}{1,000,000}$ should be the fraction of all counts that come from a gene $i$. Multiply this by the total number of reads to get the original counts.
The total number of reads may or may not be critcal, I think this might come out in the wash of the normalization, but it might mess with the disperision estimates, I'm not sure. But it is absolutely critical that
1) The TPM dataset was counted across exactly the same annotation as the count datasets
2) The gene lengths used are identical.
This could be untrue for several reasons:

Different annotation or annotation version (ENSEMBL, UCSC, REFSEQ, CCDS etc).
Different subset of transcripts used, or different way of combinng transcripts.
Different ways of calculating "effective length", used to correct for things like which parts of the gene reads can or can't map to (reads too close to ends, unmappable regions of the gene etc).

Much safer to convert the count sets to TPM, but this would preclude doing count based statisitcs (e.g. DESeq2, edgeR, voom).
However I'd also worry that if these datasets were processed at this point differently, they are going to have been processed in different ways, and that the data is basically uncomparable.

Thank you very much for the answer. Could you give any links to python (or R) packages that implement rlog? Could I apply rlog both to UMIs and TPMs? Even if not, then it might be a choice for me to do TPM for all the samples and then rlog.

rLog is comes from the deseq2 r package. It is designed counts, but as far as I'm aware it makes no assumptions that wouldn't also hold true for TPM.


text:
Check if a string is ensembl gene id or not in python
I have many scRNAseq datasets saved as hdf5 files (hierarchical data format*), and when I open them keys sometimes are confused, so I need a way of determining from string values themselves whether the array contains gene symbols or ensembl ids. What would be a regular expression in python to check for it?
*, A database format for quickly accessing large amounts of data, which is accessed using keys. The format is quicker than pickle/unpickle and can be viewed via a hdf5 viewer, a bit of info is here

Fixed it myself actually, not that hard. I found regular expression for ensembl ids here:
https://www.biostars.org/p/99142/
And used the re package in python:
import re

ensembl_regexp = 'ENS[A-Z]+[0-9]{11}'

is_ensembl = np.all([bool(re.match(ensembl_regexp, gene)) for gene in gene_symbols])
    if is_ensembl:
        gene_symbols = list(f['ensembl_ids'])

Where gene_symbols is a numpy array that I need to check whether it has ensembl ids, and if yes, switch to using gene symbols which are going to be under ensembl_ids key due to a mistake of someone who made the hdf5 dataset.

What kind of Ensembl IDs are matched by the second part of that regex? ENS[A-Z]+[0-9]{11} matches the standard Ensembl IDs, what does the rest do? It looks like it matches FlyBase IDs (CG[0-9]+), and all sorts of other things. For example, that regex would consider CG1 or ABC123bobbie valid Ensembl IDs.

Not sure, that’s what I got from the ‘Biostars’ link that is present in the answer. Maybe ask the person who answered it on ‘Biostars’ ?

Is this your question originally? Did you test this solution? The regex matches all sorts of IDs, many of which are not Ensembl. I have edited your answer so that it only matches Ensembl IDs as the question asks for.

Of course, this is still not entirely correct since it will allow something like ENSABCDEFGHIJKL12345678901 which is clearly not a valid ID.

The answer by terdon is probably good enough, but some commenters have noted it will allow something like ENSABCDEFGHIJKL12345678901.
It's helpful to reference the official documentation of Ensembl Stable Ids:

Stable IDs are created in the form ENS[species prefix][feature type prefix][a unique eleven digit number]. For example a mouse gene might be ENSMUSG###########. This means that we can immediately tell from a stable ID what kind of feature they refer to and what species they are in. The prefixes are listed on our prefix page.

The two places I look first for regexes are identifiers.org and wikidata. If these are wrong, best to report that upstream. Here are their ensembl regexes:
The regex used by the identifiers.org registry (for any ensembl ID) is quite long:
^((ENS[FPTG]\\d{11}(\\.\\d+)?)|(FB\\w{2}\\d{7})|(Y[A-Z]{2}\\d{3}[a-zA-Z](\\-[A-Z])?)|([A-Z_a-z0-9]+(\\.)?(t)?(\\d+)?([a-z])?))$

Meanwhile Wikidata has this regex for an Ensembl gene ID:
(ENS(|MUS|RNO)G\d{11})|(Y[A-P][LR]\d{3}[CW](-[A-G])?)


text:
Finding gene name from human genome using SP1 transcrition factor binding site from Postion Weight Matrix
I have a list binding site motif of SP1 transcription factor which collects from PWMScan tools. Actually this tool scan whole human genome and give the entire binding motif list. From that list of binding motif I want to align with human genome so that I can find the corresponding gene name.
If any body prescribe the suitable tools name or it will require for any programming.

Welcome to the site. Could you please tell us if you know the position of the genome of those binding motifs? Or in which format do you have that list of binding motifs?

Can you post the first few lines of your list so that we know what format it is in? Do you want to find gene promoters that contain this motif?

yes ....I am giving you two line of the format.. chr1 875622 875633 GCCCCTCCCCC 1676 - MA0079.3 P-value=8.34e-08
chr1 877443 877454 GCCCCTCCCCC 1676 - MA0079.3 P-value=8.34e-08

yes ....I am giving you one line of the format.. chr1 875622 875633 GCCCCTCCCCC 1676 - MA0079.3 P-value=8.34e-08                                                                       this is for chr 1      start site end site     the sequence ....then P value.........     Yes I find the gene promoter that contain this motif.     If you kindly recomend  me  tools name then it is very help ful for me.I found this binding motif usning PWMScan tool(https://ccg.vital-it.ch/pwmtools/pwmscan.php)

In BED format  the list of binding motifs@LIopis


text:
Finding URL for UCSC data hub
I am trying to build a UCSC track hub (a private one).
I have all my hub files on a university server in a simple folder called /hub. I am using Ubuntu 16.04 and the server is also Ubuntu (16.04.3 LTS (GNU/Linux 4.4.0-98-generic x86_64).
I have tried using hubCheck and also simply upload the hub as instructed in the UCSC instructions.
The URL I have used to indicate the location of my hub is:
http://pcxxx.xx.xx.ac.uk/home/username/projects/project-name/hub/hub.txt

However, in both cases, I get the following error:
Couldn't open http://webserver-id/home/username/projects/project-name/hub/hub.txt
TCP non-blocking connect() to webserver timed-out in select() after 10000    milliseconds - Cancelling!

I am not very knowledgeable in servers and protocols (http, ftp etc). But am I right to think that this error means that my hub is somehow not internet accessible? I'm struggling to understand what this error means. When I log onto the server, I automatically arrive in my home folder: /home/username. My hub is in /projects/project-name/hub.
Is there a way I can check if my folder location is indeed internet accessible?

Please [edit] your question and explain how you've set up your server in more detail. What operating system are you using? What server software? When do you get this error? What instructions did you follow (your link doesn't show any instructions; tell us what you did). The more details you give us, the likelier it is we'll be able to help you.

Apologies. I have just added more details.

Thanks, now we know the operating system, but the other questions are also relevant. Please edit and clarify. Also tell us whether you actually use http://webserver-id/home/username/projects/project-name/hub/hub.txt as the URL. I assume you use an actual URL, right? We need to know what kind of web server you are using (apache?), how you've configured it and all the other things I asked in my previous comment.

So what is webserver-id? Are you using the actual string webserver-id? not the ID of your web server? If so, you need to take a step back and learn about how to set up a webserver. I suggest you read https://help.ubuntu.com/community/ApacheMySQLPHP or other tutorials (search for LAMP). If I misunderstand, please edit and clarify what works, how you set it up and what should work but doesn't.

Just added the URL I actually used. As I said, I have no real idea if this is how you contruct a URL. As far as I know, the server is using Canonical landscape.

No, I just didn't put the real webserver ID. The webserver is called something like: pcxxx.xx.xx.ac.uk.

Check your webserver functions from the top down.  You can do this from any WEB browser, even one that is on the webserver system itself.  For example, does this link function:
  http://webserver-id/
and provide you with a home page for that WEB service ?  If not, then there
is not apache WEB server on that system, or is mis-configured.
Next, depending upon how that Apache is configured, try this URL
  http://webserver-id/~username/
which would be a normal way to get to user directories via Apache.
This would actually be looking for a special directory in your HOME
directory, perhaps something like $HOME/public_html/
You would need to place something in public_html/ to access, for example
a tiny WEB page:
<title>my home directory</title>

hello world

If that works, then you need to place your hub files in this
public_html/ directory.


text:
How to automate NCBI genome download
I need to download all the completely assembled cyanobacterial genome's GenBank file(.gbff) from NCBI(RefSeq or INSDC ftp data).
For this I think, the steps are:

Need to find the completely assembled genomes.
find the GenBank file URL based on the taxonomic name.
download the GenBank file(.gbff file)

Is their any way to do this by using any python module or any other idea??

Biopython has tools that could help you for steps 2 and 3, but I'm not sure how to get the information of your step 1.

can you please tell the tool name for step 2.

@Arijit Searching for "Download biopython Genebank" lead you here

Outline of solution:

get this file: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt
filter for lines having "Complete Genome" in column 12
filter for lines having a taxid (column 6) that is a descendant of taxon id 1117 (phylum Cyanobacteria), you can use nodes.dmp from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
The file download URL is the concatenation of column 20 and the last field of column 20, after separating it at '/', followed by appending '_genomic.gbff.gz'.  For example: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/007/365/GCF_000007365.1_ASM736v1/GCF_000007365.1_ASM736v1_genomic.gbff.gz

Edit:
Since I will need something like that soon, I made a Perl script for downloading genomes.
For Cyanobacteria, you would do ./download-refseq-genomes.pl 1117

As an alternative, you can use the NCBI Datasets CLI, in a two step process to generate a dehydrated dataset for later hydration:
# Download metadata for all Cyanobacteria
# And pull out accession for 'Complete' genomes
./datasets assembly-descriptors taxon Cyanobacteria --limit ALL \
     | jq -r '.assemblies[].assembly | "\(.assembly_accession) \(.assembly_level)"' \
     | grep Complete \
     | cut -f 1 -d ' ' > accs.txt

# Download a dehydrated dataset for the provided set of accessions
./datasets download assembly -i accs.txt -f 1117_complete.zip \
     --dehydrated \
     --exclude-protein \
     --exclude-rna \
     --exclude-seq \
     --exclude-gff3 \
     --include-gbff

You may be interested in some of the other files that are excluded.
Once you have the dehydrated dataset, you can rehydrate it when you are ready to analyze the data with:
unzip 1117_complete.zip -d 1117_complete
./datasets rehydrate -f 1117_complete

Use the NCBI Datasets command line tool to download genomes by taxon (NCBI Taxonomy ID, scientific or common name at any tax rank).
datasets download genome taxon cyanobacterial \
  --assembly-source refseq --assembly-level complete_genome \
  --exclude-genomic-cds --exclude-gff3 --exclude-protein \
  --exclude-rna --exclude-seq --include-gbff \
  --filename genomes.zip


text:
Running SnakeMake on cluster
I do not understand how to specify correctly the parameters on a SLURM cluster for snakemake to use them. I tried submitting the following SLURM file, but it does not work this way and the number of cores used is only 1, not 20:
#!/bin/bash

#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=12:00:00
#SBATCH --output=snakemake%A.out
#SBATCH --error=snakemake%A.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

snakemake

Then, I tried following snakemake tutorial. And created cluster.json based on the SLURM parameters that I need:
{
    "__default__" :
    {
        "A" : "overall",
        "time" : "24:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 20,
        "p" : "standard",
        "mem": 120000,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

And ran snakemake inside a newly created snakemake.sh script:
#!/bin/bash

snakemake -j 999 --cluster-config cluster.json --cluster "sbatch -A
{cluster.A} -p {cluster.p} \
-t {cluster.time} --output {cluster.output} --error {cluster.error} --
nodes {cluster.nodes} \
--ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem
{cluster.mem}"

And it is giving me an error now:

sbatch: error: Unable to open file
  /bin/sh: line 1: -t: command not found
  Error submitting jobscript (exit code 127):

I am now completely lost at what I shoould actually do. I would prefer plain regular .slurm file submission, but how to make snakemake use them? Any suggestions would be greatly appreciated.
I removed \ - line separators in the snakemake.sh script:
#!/bin/bash

snakemake -j 10 --cluster-config cluster.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --nodes {cluster.nodes} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

And it started to run. It is not convenient though for me. I would rather prefer submitting just one job using .slurm file passing all of the parameters from #SBATCH. Is it possible?

What does the %A mean?

username who is in charge of the cluster account computation resources, their payment. You can be this user, but it can also be that it is the other person, depending on who is paying for the cluster.

Can you use {rule} and {wildcards} in the output and error names in the slurm.json to name the output and error files of a job? so instead of "output": "snakemake%A.out", have ""output": "snakemake{rule}_{wilcards}.out"?

good question. idk, but i would also want to learn that since it simplifies things

You can omit --nodes, you need the following:
#!/bin/bash
#SBATCH --ntasks-per-node=1
#SBATCH -c threads_from_snakemake
#SBATCH -p some_partition

Commands go here

For slurm, you may want to modify my SlurmEasy script. We use that with snakemake all the time, usually of the form:
snakemake -c "SlurmEasy -t {threads} -n {rule} --mem-per-cpu {cluster.memory}" ...

You'll also want units with your memory requests, like 12G. Also, in general it's best not to submit snakemake itself to the cluster, rather run it on an interactive node and have it submit jobs on your behalf (if the admins complain point out that snakemake is barely using any resources on the head node).

While I agree with not submitting snakemake to the cluster, it's not a great solution to make your pilot program run interactively - automation via batch script should not be dependent on an open shell.
How about allocating a whole node with many cores and working within the GNU parallel abilities of snakemake.


text:
bioawk removed part of FASTQ header
I used
bioawk -cfastx 'length($seq) > 1 {print "@"$name"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

in order to keep a particular read length, but this command shortened the header  from
@A00199:161:HF3JLDMXX:1:1101:5882:1063 2:N:0:CTGTAGCC+GAACAGGC

to
@A00199:161:HF3JLDMXX:1:1101:5882:1063

How is it possible to fix the bioawk command to remain the full command?

bioawk -cfastx 'length($seq) > 1 {print "@"$name $comment"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

Note the addition of $comment, which is what follows the name on the first line of each entry.

Space was missing and I had to change it to {print "@"$name " " $comment"\n"$seq"\n+\n"$qual}


text:
Counting the characters in a row using python
I am using python 3 and biopython 1.72. I have been using a protein sequence file in clustal format.
CLUSTAL 2.1 multiple sequence alignment

Human           ------------------------------------------------------------
Chimpanzee      ------------------------------------------------------------
Dog             ------------------------------------------------------------
Mouse           ------------------------------------------------------------
Xenopus         ------------------------------------------------------------
Amphioxus       MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA

Human           ------------------------------------------------------------
Chimpanzee      ------------------------------------------------------------
Dog             ------------------------------------------------------------
Mouse           ------------------------------------------------------------
Xenopus         ------------------------------------------------------------
Amphioxus       HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK

Human           ------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS
Chimpanzee      ------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS
Dog             ------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS
Mouse           ------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS
Xenopus         ------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST
Amphioxus       INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD
                                  *:***::  :.  :::     *** **: *        * .

Human           NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI
Chimpanzee      NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI
Dog             DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL
Mouse           DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL
Xenopus         EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ
Amphioxus       DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH
                :  * *:** :*.* **::.* . **:***** ::* .   .               *

I need to iterate over the rows of my file so that I could count the occurrence of characters in each row.
I tried using the code:
file_var = input("File: ")
count = 0
with open (file_var , 'r') as fo:
    for i in fo:
        count += 1
        print (i)
    print (count)

But this results just counting the total number of rows in a file(count variable). If I iterate over the file using biopython, it does'nt work also.
for seq_record in SeqIO.parse("sample.aln","clustal"):
    print (seq_record.id)
    print (repr(seq_record.seq))

seq_record.id results the names of the sequences used in the file.
repr(seq_record.seq) results the sequences of each specie alongwith description.
While I need to

get the rows only with amino acids (exclusing the symbols '*,.,:' line) and nothing else from the row and for each row I need to count the number of characters (exclusing '-' dashes)
the number of characters of that row be placed at the end of each of the row.

for example: for above example file, the counter should work like;
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42

you can count for the character using 'length' function and before the counting character, use regular expression to remove the other character from the string.

yes it does, but it won't do for each row separately. Thankyou for the response. :)

I think biopython is not the right tool if you want the counts row-by-row: biopython treats the sequence of a taxon as one unit, regardless of the original layout in the alignment file.

Here's one way with pure Python:
#!/usr/bin/env python3
import sys
import re

with open (sys.argv[1] , 'r') as fo:
    for line in fo:
        fields = line.split()
        if len(fields)==2:
            print(fields[1],len(fields[1].replace('-','')))

Running this on your file gives:
$ foo.py file.aln
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS 42
------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS 42
------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS 42
------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST 42
INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD 60
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL 45
EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ 45
DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH 60

Some other options:

awk
$ awk 'NF==2{l=$2;gsub(/-/,"",$2); print l,length($2)}' file.aln
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS 42
------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS 42
------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS 42
------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST 42
INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD 60
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL 45
EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ 45
DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH 60

Perl
$ perl -lane '$l=$F[1];$l=~s/-//g; $#F==1 && print "$F[1] ", length($l)' file.aln
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS 42
------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS 42
------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS 42
------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST 42
INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD 60
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL 45
EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ 45
DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH 60

Well ,Thankyou for this much help, it really solved my problem. Would this also work for a phylip file format?? which is ; [phylip format example] (http://scikit-bio.org/docs/0.2.3/generated/skbio.io.phylip.html)

@SidraYounas I'm afraid not. This approach assumes that the lines of interest have two fields. Phylip has more. You might want to ask a new question about that.

,alright. But how come I write up these countings into the file and not just printing them out.?? As I'm working with windows and not linux system

this code does not write the character numbers in my file. Do I need to use file handler to access writing permission to file??

@SidraYounas yes, that's just normal python. Next time, please make sure to include all your requirements in your question so we don't need to do this sort of back and forth in the comments.

I tried like this:with open ("sample_clstl_seq.aln" , 'r') as fo:     for line in fo:         fields = line.split()         if len(fields)==2:             fo.write('{0}{1}'.format(fields[1],len(fields[1].replace('-','')))) But errors out that it is not writeable.

I understand, I will be more concerned about this matter.

@SidraYounas I would not advise modifying the file directly, unless you have some good reason to. It is simpler and less risky to generate a new file and replace the original once you are sure that the result is correct.


text:
Output file names are not known beforehand
I have a folder where the hdf5 files inputs are (matrices of scRNAseq raw UMI counts). There is a python script that loads them, splits them into multiple matrices, normalizes them, and saves as csv files. When the name pattern of the output files is known, things are easy, I know how to do that, but when one input (hdf5 file) can produce various number of outputs (csv files) I am not sure how to approach the issue. Currently my Snakemake file relevant part looks the following way:
hdf5_files,=glob_wildcards('/main/{hdf5_file}.hdf5')

rule all:
  input:
--->  expand('data_files/normalized/{hdf5_file}.csv', hdf5_file=hdf5_files)
    ...

rule normalization:
  input:
     script = 'python/normalize.py',
     path = '/main/{hdf5_file}.hdf5'
  output:
--->  path = 'data_files/normalized/{hdf5_file}.csv'
  shell:
     "python {input.script} -data {input.path} -path {output.path}"

Above I highlighted the 2 lines that I am not sure now how I should write. Any suggestions would be greatly appreciated.

Have you considered using flag files? Or do you need the names for downstream processing too and can't use glob or specify a directory at that point?

The rule normalization processes, say, 2 hdf5 files and outputs 10 csv ones into data_files/normalized folder. Then, I want to use glob_wildcards to read in these 10 csv files and move on.  Further rules are, of course, fully dependent on the rule normalization because it generates the input csv files. I am not sure how can I use flag files, but will try. Thank you!

Ah, if you need to process each of the resulting csv files independently downstream then there are no snakemake tricks that I know of that will help (I presume glob_wildcards() is run during rule parsing, so the CSV files won't be there yet). I suspect you'll want to run normalization from a wrapper script beforehand in that case, since then globbing will be certain to work as desired.

What do you mean from a wrapper script? Outside the pipeline, you mean?

Yes, exactly. A (probably python) script that runs two Snakefiles when you call it, with appropriate options. Then you have the flexibility you need while maintaining the convenience of snakemake.

Could you give an example of names of input files and output files of normalize.py ?
It seems to me that your case could be similar to the one described about the dynamic output:
https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#dynamic-files
Here you would probably use it for the output of normalization and for the input of all (with a wildcard different from hdf5_file).
Disclaimer: I haven't had the opportunity to use this feature so far.

Sample name is {hdf5_file}_{num}.hdf5, where num is an integer.


text:
Best data repository to publish a large 'plotly' table containing all annotations on a transcriptome?
I am not sure whether this is the best SE community for this question.
I want to publish a transcriptome paper, along with interactive materials enabling readers to peruse the data behind the discussion. The R package 'plotly' enables rendering online-publishable tables that can be filtered & searched for specific entries.
I would like to publish such a table as supplementary to a paper (e.g. preprint version) however usually data repositories allow for uploading files and text.
Would anyone here know an official data repository with a format such that would enable publishing online an interactive table?
A good example of such a table is given here.
Thanks in advance.

This might be better suited to [Academia] where people from more fields could know about a place for such things.

@Llopis I will also ask there, thanks! My guess is that they'll think it is mostly technical and redirect me elsewhere. I think it is in fact a quite specific question on a situation not many have faced yet.

You may try Code Ocean. Currently, it is not that easy to use IMHO, but one of their tech leads promises to substantially improve it towards the year end.

I have just looked into Code Ocean as you suggested, and though it looks promising it is indeed not very user-friendly at the moment... Plus, looks like one needs a subscription to run interactive content.

I have seen some (1 or more) articles with interactive data displays in the F1000Research journal I don't know recall now if they had interactive tables or plots but they used plotly for sure.

Yes, I heard that F1000 encourages plotly content. I have just sent them an email. However that would be a journal option for a final publication (if I can pay their fees). Ideally, right now I'd greatly benefit from an external repository that could also be used with a preprint server like Bioarxiv. Unfortunately bioarxiv seems to only take .pdf output.

Yes, I realized later that your data seems "huge" so a publisher doesn't seem a good choice to store the data

Mind that by "huge" we mean a table is very large in numbers of rows. In reality the actual data behind the table are .txt files adding up to less than 20 Mb. I could in principle publish such table as a .pdf file, however nobody would like it and read it.

The example looks a lot like a shiny application. Your best bet then is one of the following:

Ask your IT department to setup a shiny server for you. They can use the free version.
Post your shiny app on shinyapps.io. Depending on how popular the page becomes, you might have to pay for an account.

The free version of shiny works for most cases (it's what we use), the only problem with it is that it doesn't scale (but if you care about that then you need to either pay for hosting on a pro server or buy a shiny pro license).

If you have problems scaling shiny apps you might find this post interesting.

Only now I see the issue you mention about popularity. On a free account I only get 24 mins of server availability per month. Given my scope of use, this is way too short !..


text:
How is prior.count used by edgeR's cpm
edgeR's cpm function has an argument called prior.count. Based on my understanding of the documentation, it is supposed to be adding a fixed number per sample which is proportional to the library size of said sample. Average of all numbers added to all samples would be equal to prior.count.
However, looking at actual data this does not seem to be the case. Given an imaginary data frame of
df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)

We can try to calculate log cpms by doing
logCPM = cpm(df,log = TRUE,prior.count = .5)

We can also calculate regular cpms by doing
CPM = cpm(df)

If we wanted to see, what is the number that is added to the CPMs before getting logged, we could do
difference = 2^logCPM - CPM

To make it look pretty lets use tibble
tibble::as.tibble(difference)

# A tibble: 4 x 3
                     a                    c                    b
                 <dbl>                <dbl>                <dbl>
1  25641                25641                25641
2  12821                12821                12821
3 -    0.0000000000582 -    0.0000000000582 -    0.0000000000582
4  38462                38462                38462

Here we see that the number that is added has no relation to the library size of the sample. What I want to learn is, how is this number that is added to each cpm is calculated based on prior.count.
I have been digging through the code but that part goes into C++ territory which hinders my understanding of what is going on

I get different results in the difference table. Which version are you using ? (I am in R 3.5.1 and edgeR 3.22.2)

The prior count ends up getting scaled by the ratio of a library size to the average library size and then multiplied by 2 before getting added to each library size (I'm sure there's a good reason for that, but I don't know what it is). Using the example df data frame in your post, let's walk step-by-step through what cpm() is doing:
df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)
prior.count = 0.5
# First, we need to calculate a library size
lib.size = colSums(df)
# Calculate the average library size and the adjusted priors
ave.lib = mean(lib.size)
adjusted.prior = prior.count * lib.size / ave.lib
# Update the library sizes
adjusted.lib.size = lib.size + 2*adjusted.prior
# Now we can compute the CPM
customCPM = t((log(t(df) + adjusted.prior) - log(adjusted.lib.size) + log(1000000))/log(2))

The matrix transposition looks really ugly in the last step, but it's needed so R adds things like adjusted.prior to each row of df.
Note that the prior count is only used if you compute log(CPM), it's unused otherwise.

How do you know that edgeR is scaling that way to obtain the adjusted.prior?  By the way sweep  could be used to avoid having to transpose the matrix

It's stated in the source code. I know I can use sweep(), but t() takes fewer characters.

Thanks a lot. By the way was I wrong about the code being within the C++ part? If that R code above was part of the package, I have missed it

You didn't miss it, I translated the C++ code to more accessible R. The C++ code is pretty tough to follow unless you're pretty comfortable with it.

thanks this is an awesome answer and helps me understand prior.count


text:
How to map short sequences to long reads, recovering all multiply-mapped high-quality matches
The dilemma:
I have around one million short sequences (21 bp to several 100s of basepairs) for which I need to identify all occurrences of in 20-30x coverage noisy long reads (both pacbio and ONT). All of the short sequences and long reads are derived from the same individual, and there are enough short sequences given the organism's genome size that each long read should have at least one or multiple short sequences.
An example:
Given three short sequences, AA & BBBB & CCCCCC , how to identify all occurrences of these sequences in long, uncorrected reads.
Sequences to identify in long reads:
  - AA
  - BBBB
  - CCCCCC

Long Reads:
  - 1: ------AA-----------------------CCCCCC------
  - 2: ------AA--------BBBB-----------------------
  - 3: ---BBBB-----------CCCC--
  - 4: ------------------AA-----

In this example, the output sam file will contain high MAPQ hits for short sequence AA to reads (1, 2, 3), for short sequence BBBB to reads (2, 3). For short sequence CCCCCC there will be a high MAPQ hit to read 1 and a lower MAPQ hit to read 3
Limitations:

Due to the requirements of this project, I cannot perform a de novo assembly then map the short sequences to the assembly output. It is necessary to directly map the short sequences to the long reads.
I cannot correct the long reads before I map/align the short sequences to them.
The alignment/mapping should sensitively detect the short sequences in the long reads.
The short sequences may be anywhere from 1Mbp to 30Mbp in total bases.
The long reads may be anywhere from 3Gbp to 90Gbp in total bases. So the mapper/alignment technique is hopefully fast. Sensitivity is more important though.

Research so far:

bwa mem does not output multiply-mapped in a predictable manner that fits this use case https://www.biostars.org/p/304614/
blat seems to output many hits for a single short sequence, but there are no out-of-the-box parameters for mapping short reads to long reads.
Maybe multiple sequence alignment is the most sensitive, but wouldn't that entail running (no. of short seqs) * (no. of long reads) alignments?

Don't map short reads to noisy long reads. This won't work. As you already have decent coverage, assemble first (or correct long reads by themselves with canu) and then map. It will be faster and give you much better result.

Hi user17818 - unfortunately that isn't an option in this scenario based on the requirements for this project. I'll update my question to better clarify that de novo assembly, then mapping the short sequences is not an option.

This is probably an intermediate problem of a larger project. If so, you need to change the strategy. Mapping a 21bp oligo to 90Gb noisy reads will lead to many false positives and negatives. This is a theoretical limitation that no direct mapping algorithms can get around.

You're right - there will be false negatives and positives. My goal is to minimize them by finding the best mapping strategy to get multiple mappings of each query. The size range of the illumina-derived "short" query sequences vary from 21 to hundreds/thousands of bp. The chances of correctly mapping the short queries to a long read will increase as the lengths of the illumina-derived queries increase. The kmer space of even a 31-mer is 4.61x10^18, or around 9 orders of magnitude larger than the kmer space of a human genome. A few indels in the long reads should still be uniquely mappable.

"Uniquely mapped" doesn't mean the mapping is correct. This is particularly true when you look for all partial mappings. For >100bp sequences, you may use minimap2, taking Illumina sequences as the reference. It won't work well with shorter sequences, though.

LAST has given the best results for me when I've tried to do this, although I agree with @user172818 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see here), as well as abundant base calling differences present in single-molecule sequencing.
Minimising error is not necessarily going to be the best option, and concentrating only on unique hits will miss a lot of signal. There are frequently multiple identically plausible positions, even at zero error, and the long reads will have their own associated errors.
I also find the limitations a little odd. Canu can do read-level correction on long reads based on overlapping reads, and if assembly is possible, then nanopolish can be used to correct nanopore reads for systematic base-calling error introduced by the methylation of unamplified template.

Thanks for your comments, @gringer. The limitations in the problem are the limitations. The best I can do in this case is minimize false positive hits. I'll see what I can do with LAST!


text:
Article GEO gives counts that are not integers. Should UMI counts be integers?
I thought that UMI counts are always integers, but when I opened several datasets provided by GEO I got confused because in some articles they provide counts that are supposedly non-normalized raw UMI count values, but they are not integers. For instance: GSM2043231 or GSE77067
Could anyone clarify, what these people mean by abundance measurements?

The data processing section doesn't mention using UMIs at all, so it's likely there were none (whether they'll be present depends on how the libraries are made, with lower throughput methods being less likely to include them). The values are stated as being expected counts from RSEM, which is why they're not always integers.


text:
Get genomic coordinates using GenomicFeatures by HGNC gene names
I want to get coordinates of human genes from my list (consisting of hgnc genes id) using GenomicFeatures and TxDb.Hsapiens.UCSC.hg19.knownGene R packages from Bioconductor.
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb=(TxDb.Hsapiens.UCSC.hg19.knownGene)

my_genes = c("INO80","NASP","INO80D","SMARCA1")

select(txdb, keys = my_genes,
   columns=c("TXCHROM","TXSTART","TXEND","TXSTRAND"),
   keytype="GENEID")

However, it doesn't' work because txdb doesn't take hgnc identifiers; how can it be solved? I couldn't find any appropriate keytype that will support hgnc and not sure how to match hgnc id I have and GENEID from txdb.

Answer from @devon-ryan, converted from comment:
Try a UCSC gene ID instead. For those genes they'd be: uc001zni.3 (INO80), uc002vaz.4 (INO80D), uc004eup.4 or uc004eun.4 (SMARCA1) and uc001coj.2, uc010olq.2, uc010olr.2, uc001coi.2 or uc001col.2 (NASP). You might have to omit the versions (the .4 part, for example).