text:
MaCconkey Agar colour change in non-lactose fermentors
On Macconkey agar, lactose fermenting bacteria turn the agar bright red. Why does non-lactose fermenting bacteria turn the MaCconkey agar yellow?. Does it have to do with the alkalinity of waste products?

I'm voting to close this question as off-topic because it belongs on Stack Exchange Biology. There is nothing wrong with the question, just wrong site.

McConkey agar is already red. It's trick is that it has two carbon sources. The bacteria that cannot use the more energetically useful lactose use the amino acids (peptone) generating ammonia which changes the pH and the crystal violet goes yellow. The ones that can use lactose acidify the agar a bit so go a bit more red —but not much.
This being the SE for bioinformatics, it should be said that if the genome of the organism in question is sequenced and deposited in NCBI you can find if the organism possesses the enzyme that allows the hydrolysis of lactose (beta-galactosidase) by doing a protein BLAST. The lacZ gene encoded protein in E. coli is NP_414878 in NCBI. This goes in the accession number input box. The species in question goes in the organism box. If you get something with good coverage (80%) with more than 20% homology the organism is likely lac+.


text:
miRNA dataset for prediction experiments
I am from the field of Computer Science and I want to perform some experiments related to miRNA and mRNA prediction.I have been searching for datasets in this area and found one called mirWalk available at: http://mirwalk.umm.uni-heidelberg.de/
when I enter, for example, the Entrez ID 10608 it appears to me the following table (I have only put the first record):

I have two questions:
a) Can I use this dataset to try to make predictions about the interaction of this gene with, for example, miRNA hsa-let-7c-5p? For what I see this interaction appears in the table, but it seems it has not been validated.
b) What are the meaning of the different fields, there is no information about it in the website of mirWalk?
Thanks for your help and consider that I am a newbie in this bioinformatics field.

source: Stack Exchange

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

sort: votes


text:
Clusters in scRNA of Glioblastoma
There is a scRNA dataset on Glioblastoma from 10xgenomics.
Unfortunately, there seems to be no research paper related to that data.
Question 0 : May be there is paper  ? How to find it ? I am not familiar with 10xgenomics - may be there is some standard way to find papers related to their datasets ? May be their data also placed in GEO database, then what is GSE-id ?
Any way, we can see several clusters in the data (see figure above or switch to tab 'analysis' at the link above).
Question 1: What can be biological interpretation of these clusters ? I would expect cancer cells to be more homogeneous, but here we see several well-separated clusters, probably something like that has been already discussed in the literature ?
PS
That dataset has been placed on kaggle https://www.kaggle.com/reminho/human-glioblastoma-dataset
by enthusiast
and one can find several notebooks in R and Python to play with it. Clusters are also present there.
In general there are many datasets/notebooks/people related to bioinformatics on kaggle, it is quite convenient platform for everything related to data science - data sharing, data processing on their servers for free (like colab), commenting on the works of others and so on.
It seems it would be useful to create bioinformatics "community" on kaggle (there are already 12 other "communities") if you find that idea - useful please share your comments/make upvote/show yourself at:
https://www.kaggle.com/general/203136

I actually do not think cancers are homogeneous. There can be cancer cells, normal cells, various immune cells in cancer tissue. So it is not surprising to see these many clusters, it may be helpful to try to find markers for different clusters and identify the cell type.

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


text:
Overlap in Paired-end Reads for Sequencing?
Regarding Sequencing:
Do/can the paired-end reads have overlaps?

Yes, if the insert size is smaller than the read sizes, this would happen. For some applications, for example SNP detection for molecular diagnostics purposes, this approach is used.


text:
Location of Reads in Sequencing
I have questions regarding sequencing:
Are paired-end reads from different stands?
What about single reads, are they coming from both strands?

Mates in a pair are giving signal for a single strand, though this is only apparent in bisulfite-treated samples. The reason for this is that you are sequencing the ends of fragments of material and it's random which strand of this (or both in different clusters) will get sequenced. Different pairs of reads will then randomly sequence one strand or the other.
Single-end reads are the same as paired-end reads, the second read in the pair is just not made.

Now, I understood from you that the mates in paired-end reads are coming from two different strands, but based on the figure in https://www.biostars.org/p/95803/ and also some other sources, I implied that they are coming from two different stands as they are read from two different directions. Could you please help me in understanding this?


text:
Warning in fastqc
I am checking the quality control of my sequences using the Fastqc tool. For some steps, like "Per base GC content", etc., I received a warning. So, I was wondering whether I Should also take care of warning messages in the Fastqc tool or I should only care about errors!
I would appreciate it if someone could help me with this, please.

Whether the warning is heeded or not it all depends on the kind of data you are working on and the type of analysis you should be doing. If you could tell us a little bit about what you are working on and what kind of analysis you want to do that could help us

Generally yes, but it depends the type of study. Remember that Fastqc was developed for bacterial genomes, mainly, so an unusually shaped distribution in that particular dataset could indicate a contaminated library. But, for example, in an amplicon 16S amplicon dataset it's expected that other modules fail (Duplicate Sequences, Overrepresented Sequences, etc.) not only GC.

Maybe, it is good to mention that I am working on the genomic data of E. Coli. So, in this case, you believe that the warnings are also important, right?

According to the documentation "A warning is raised if the sum of the deviations from the normal distribution represents more than 15% of the reads... usually indicate a problem with the library". You need to check in the downstream analysis if you have adapters, or other kind of contamination.


text:
Can ortholog genes have different names among different species?
To be specific, we have different lineages of E. coli. I wondered if the ortholog genes among them should all have the same name in all lineages or can they also have different names?

Gene naming is pretty subjective and even a single gene can have many names. If you're lucky orthologs will be named the same in different lineages, but you can't rely on it. Better to use gene ids (e.g. from Ensembl) and an ortholog db to map the orthologs. However, if you're looking at different strains, this is more tricky since most db don't go below the species level

Thanks for your answer. I was wondering if two genes belong to two different isolates of E. coli (i.e., they are still the same species but different isolates), should we call them ortholog or paralog

Gene IDs are not perfect either. Those can change as well, even between different versions of the annotation for the same species.


text:
Evolutionary speaking, what are the constraints on having new mutations?
Evolutionary speaking, some mutations lead to better fitness of an organism and its adaptation to the environment changes, but I was wondering if there are some constraints on having new mutations. For example, if a gene is needed for survival, I expect to see it more resistant to having new mutations. So, my question is what are these constraints, especially for bacteria?

I think the is a very broad question, do you have a more specific problem in your mind? Or example you are puzzled by?

This is an undergraduate/graduate lecture type question and the answer varies: traditional Darwinists believed a reasonable proportion of mutations would result in beneficial adaptation, neutral theory however would state the vast majority of mutations have no effect on phenotype and only an absolute minority result in beneficial adaptation. It depends on which bacteria are the focus of investigation and the selective constraint (I would assume its antibiotics).
Once the bacterial species and selective constraint are described its a straight-forward empirical investigation using molecular mathematics because the mutations driving antibiotic resistance are all known.


text:
Per Base Sequence Content in fastqc
I have a question regarding "Per Base Sequence Content" plot for "fastqc":
In the fastqc documentation, it is written: "In a random library you would expect that there would be little to no difference between the different bases of a sequence run, so the lines in this plot should run parallel with each other."
But I don't understand why different bases in a read should follow the same pattern of allele frequency ("A/T/C/G"). I mean they are different positions in the genome and it is normal that each position has different allele.
I would appreciate it if someone could help me, please.

Heads up: This question has been cross-posted on biostars

The positions you see in the x-axis of this plot are positions on your reads, not on the reference genome. Since each read comes from a random position on the genome, the frequency of A/T/C/G on each base position should reflect the base composition of the entire genome. Therefore, the percentage of each base should be approximately constant across the read.
Hope that helps!


text:
Since every human has a different DNA (different combinations of C, G, A, T) what does it mean to have the genome done?
I'm confused about the difference between genome and DNA. Is it correct to say that the same type of bacteria has the same DNA? But my understanding is that it is not correct to say that the same type of human has the same DNA, since every human has a different DNA. What am I missing here?

I'm voting to close this question as off-topic because it is a basic biology question and is not really bioinformatics. Biology SE would be a better place to ask this.

Closing is never a good idea, unless it is something outrageously wrong or off-topic. Not the case here.

As the OP, you'd surely have a bias against closing your own question.

@RamRS, It's clear to me that the OP is asking what people mean when they say "the honey bee genome," "the horse genome," or "the human genome," given that most individual honey bees, horses, and humans differ genetically. Although I'm neither a biologist or a bioinformaticist, I have the impression that genome projects were deeply intertwined with the rise of bioinformatics as a field, so this looks like a bioinformatics question to me. It may be a basic one, but to me it feels basic in the sense of being fundamental.

@user253751 The current top rated answer reaches a different conclusion.

I work in genomics but am not a geneticist, so I also had this question for many years. The answers below are correct -- the human genome is a reference point from which we discover variation in others -- but it underlies how confusing genomics jargon can be. Another example is "variant calling". Why "calling"? I was told it was the meaning of "to call" as in "to call a race" and even then it took me a while to understand! Jargon is coined by metaphor, so in my opinion, genomicists must have quite an "interesting" view of the world!

All humans have some differences in their DNA, but there's far more that is shared. On average the difference between humans is only about one thousandth of their full DNA, which means we're about 99.9% the same. These differences aren't distributed fully randomly, but are often because of specific gene alternatives. (Random mutations do occur, but they are also often fatal, so the random mutations we see of living adults are much more restrictive than all the random mutations that occur.)
To identify the human genome is to study many people's DNA and to label the parts that are shared between everyone, the parts with two or three variations, and the parts with even more variations. Even though every person has different DNA, we can still say they fit the pattern, just like every T-shirt may be unique but they all fit the T-shirt pattern and not the trousers pattern.

DNA is only the substance that makes up the genome.
DNA can come with many forms. For example, the fragments you get by PCR in the lab are DNA. A bacterium also has DNA molecules in it. You may also find DNA fragments in many other places, like in our blood, or at a crime scene.
Genome, however, is a very specific term. It means all the heritable substances in an individual or an organism. Most times, when we talk about the genome in eukaryotes, we refer to the nucleus genome, that is, the DNA in the nucleus.
Note that genomes do not have to be made of DNA. RNA viruses, for example, only have RNA as their genomes. Yet we still refer to that as a genome, because that is the heritable substances for those viruses.

Might have misunderstood you - but "we" have sampled and "decoded" the DNA of a few humans (often volunteers who surrendered their body to science after their death) - but this gives of course only the genetic code of those few individuals.
When scientists say we decoded the genome they mean they deciphered which parts of the DNA "do" specific things - the actual "genes" or gene-sequences - usually done when looking at the DNA molecule packed up into the chromosomes to give easier "orientation".
And by my understanding the first task was accomplished - supported by use of supercomputers, but the deciphering of the gene-sequences still has a long way to go for scientists to state "complete." But there is growing success - And as usual scientists will not stop until complete.

Why do you think we have only/mostly sampled the DNA of the deceased? You can take a swab of your DNA (non-harmful, not even painful) and send it to various labs for identification of inherited diseases, genealogy (ancestry), etc.

I wrote about those as I wanted kinda honor those who donated their bodies - but you are right .. of course the DNA samples can come from living humans too since literally a tiny drop of blood would be enough @CJDennis

That makes you a good human being, but it's not good to have emotion in answers. Answers are best when they stick to the facts.

"To have the genome done" means that we have sequenced all the DNA in particular reference individuals. You are correct that every human has different DNA, so there is no such thing as a single "human genome", as (almost) everyone's is unique. But the vast majority of everyone's DNA is identical - that's what makes us human. Having the full sequence of a few "normal" individuals who don't have glaring genetic defects is extremely useful for learning about the genome of all individuals, since so much of it is the same from person to person.
I've used the term "reference" a few times, since we use a few individuals' genomes as an expectation of what a human genome looks like. The process of sequencing DNA involves breaking it into many bits, and then recombining them using a computer. You can put the pieces together without the reference, but it is much more difficult and time consuming. Even if there are differences from the reference, having an expectation of what >99% of the genome should be makes the problem much simpler. In a way, it's like the difference between doing a jigsaw puzzle without knowing what it's supposed to look like, and using the picture on the box as a reference. Even if your exact jigsaw puzzle has a few different pieces, the box reference is  mostly the same and immensely useful.

When people say "the genome" for a species, they mean a reference genome. To quote wikipedia,

A reference genome (also known as a reference assembly) is a digital nucleic acid sequence database, assembled by scientists as a representative example of a species' set of genes. As they are often assembled from the sequencing of DNA from a number of donors, reference genomes do not accurately represent the set of genes of any single person.

The genome is updated over time as more is learned, the most recent human reference is hg38.
For individualized results, the patient is sequenced and the data is compared to the reference genome. Then, the differences are presented as a list of genetic variants. This is more efficient than looking at the entire genome for each patient, because on average greater than 99.9% of sites match the reference.

Yes you are right.  99.9% of our genome remains same, it is only the 0.1% which makes each individual unique from the rest of us. Single nucleotide polymorphisms (SNPs) contributes to the 0.1% variation in our genome. SNPs explain the variation in eye color, susceptibility towards disease conditions etc

Think of it as a special book. Everyone's book is different, but they all have the same number of chapters and each chapter covers the same topic in each individual book (sorta). The problem is we can only read a few words at a time.
When the genome is "done", it means that someone has read a bunch of these sentence fragments from one person's book and then assembled all the words in the correct order to make a complete book. Then we can get a good idea of how many chapters it has and what those chapters are about.
However, the human genome wouldn't be a novel, it would be more like a set of encyclopedias.


text:
How can I match sequences with IDs using python dictionaries?
0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

I've asked a similar question before. I have a better understanding of python now and would like to revisit this. Can someone critique my thinking?
I have an HMM output file that contains about 20 fields. There are two sets of IDs, e.g., P006672 and jhhlp_005678
In a separate file are FASTA-formatted sequences linked to these IDs.
Proposed solution:
1
   a) read in fields containing IDs, lets called them fields 8 and 6;

   b) store field values in an array or dictionary (whatever the proper data structure is for
  this). Let's called this Fields_8_and_6

a) read proteome file containing FASTA-formatted sequences
b) Grep and store as ID:sequence key:value pairs; let's call this ID_to_sequence

a) Compare Fields_8_and_6 to ID_to_sequence
b) Write ID1:sequence1 ID2:sequence2 as pairs into one file

send to ClustalO

I'm not sure if this will be helpful to you, but my initial idea would be to use the very helpful pandas module in python to parse the table, then use some sort of regex to recover the required IDs (for example, I might set a regex query to recover 'text between vertical lines' to get P06672, and 'text after =' to get jhhlp_005678). There are tutorials for pandas and regex on the Software Carpentry website.

Got it, I know I've used both biopython and regex to parse the IDs before, and I've used awk to calculate query coverage (e.g., awk '{$1 = ($17-$16)/$6}1' ). Thanks very much for this

pandas will do it, but as pointed out BioPython has an HMM package.

Biopython can index a fastafile and you can access the sequences by their ID in a dict-like manner.
All that is left for you to do is extract the IDs from your file. (more Biopython)

I didn't know they supported parsing hmmer outputs, thank you very much

Hi Pallie, I merged my drugbank file with my proteome file; I'm getting duplicated because of the drugbank file. Is there a way around this with additional keys or something of that sort?

Hi Mike, open a new question and show your code, input, and expected output so people can help :)

sounds great thank you! This is a really powerful protocol, automating it and getting it on the web would go a long way


text:
STRUCTURE interprets pop ID as locus
I'm using structure for the first time and I fail to execute a "job". For some reason, structure interprets the pop ID column as locus.
I have following .txt file

The first column are sample labels, the second column are population IDs and then there are 26 columns containing the locus information (i.e. 13 markers). Thus, in total 28 columns. It is a diploïd organism so each locus consists of 2 consecutive columns (as described in structure manual). In total there are 64 rows (so 64 samples)
I load this .txt file into structure, specifying following parameters:
-number of individuals 64
-Ploidy of data 2
-number of loci 13
-missing data value -9
-special format: data file stores data for individuals in single line
-individual ID for each individual
-putative population origin for each individual
The loaded file looks like:

This looks good. My samples have a separate column, pop ID is a separate column for possible populations an then for each locus there are 2 columns. Not all loci are visible on screenshot, but there are 13.
I want to calculate Delta K and therefore I want to run a job. However when I run it, I get following error:

I get this error for each row/individual. It appears that structure interprets the pop ID values ("Noord","Zuid", "onbekend") as 14th locus, but I have never specified a 14th locus.  I fail to understand why this happens.
When I run the program without pop ID column, it works.
Could anyone please advise me on how to solve this error?

Hmmmm .... you don't actually need pop ID to work for Structure, providing their are no logical errors (you get a wrong answer). The bit that must work is the a priori number of populations in the data set. You are then going to manually (or via a script) feed a range of a priori values into Structure, run each and plot the likelihood. Where the likelihood plateaus thats the number of populations contained within your data. This is the key output of Structure.
Personally, I would perform a Structure run on the test data and remove the pop ID column and check the answers (likelihood -lnL) are identical. Thereafter use an approach such as $F_{ST}$ to interrogate the exact members of each population delineation.

I think your test will work .... The only thing I can think of is the header comprises a space, This format just needs checking against the test data.

Hi Michael, thank you for your answer. I checked the text file in notepad ++ showing all characters, but nothing seems wrong. I would also like to make a triangle plot in structure, isn't it therefore required to have pop IDs? Anyway, it should work even with pop IDs included. My plan for now is to run the structure analyses (from K=1 to K=9) and then analyse the data with structure harvester to determine the number of populations (delta K)

I agree that pop ID would be useful, and if membership of a population is not intuitive (e.g. geographic foci) - its imporant, but debugging this remotely is tricky. Your strategy K1 to K9 is absolutely correct. Structure is rarely used in isolation however, hence population membership can be explored notably by Fst using a variety of methods.

Yes, my next step would be to analyse the data in popgene32. Thank you for the help.


text:
local BLAST error: BLAST Database error: Error: Not a valid version 4 database
I want to use a local blastn command to BLAST a multifasta file ~700 sequences using following command: blastn -db nt -query fasta_all.fasta -num_alignments 2 -out fasta_blasted.txt but I receive error :BLAST Database error: Error: Not a valid version 4 database.
I use a local nt database which I downloaded today (11th of May 2020) by running update_blastdb.pl --decompress nt [*]. I downloaded this but I did not do any post-processing and just the "raw" download files are in the folder (see image below which shows only a part of the files)

The blast+ I downloaded and installed is ncbi-blast-2.10.0+-x64-linux.tar.gz from https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
So, both database and blastn version should be the latest. It probably has something to do with either the version of the database or the software, but as it both are the latest version, I find it strange that these would not be compatible. Any ideas what causes this error?

first I recommend you to check your downloaded data with md5 to be sure there was no problem with that step. How many files did you download?
Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called is a version corresponding to the database you are using.

Thank you for your reply. I downloaded 24 files in total, from nt.00 till nt.23. For every nt file there are following extensions available: .nhd .nhi .nhr .nin .nnd .nni .nog .nsq. For every nt file there is a nt.xx.tar.gz.md5 file available. Does that make sense?

Try to run md5sum nt.xx.tar.gz.md5 for each xx in 00 to 23. Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called from your latest version folder.

I think this last issue you adres could be the problem. When I run blastn -version it returns Package: blast 2.6.0, build Jan 15 2017 17:12:27 so not 2.10 which I downloaded. Once I also downloaded blast with conda so perhaps that is the problem.

Thomas, indeed this was the problem, very good suggestion. I removed my BLAST through conda remove blast to remove the old version. I followed the exact steps of https://www.ncbi.nlm.nih.gov/books/NBK52640/ but this did not lead to a correct installation. The command blastn was not found even after right configuration (check by echo $PATH). Now I recall that I install blast also through conda because the description of NCBI did not work. I now managed to install blast 2.10 by following these steps: https://ngsprimer.wordpress.com/. I receive no error now. Thanks!

Happy to hear it solved the problem! good luck!

I run into this problem and indeed it was mismatch of database and blast versions. So other don't have to fish it out of the comment section to, I added it to the answer (you as much your words as possible). Feel free to adjust it, but that's the solution that should be accepted I think.


text:
"megahit: command not found" after correct installation with miniconda
I want to use MEGAHIT (git: https://github.com/voutcn/megahit) to assemble metagenome shotgun sequencing reads (illumina).
I install the package with the command conda install -c bioconda megahit as described in on the github page which yields:
Solving environment: done

## Package Plan ##

  environment location: /home/datura-workstation/miniconda3

  added / updated specs:
    - megahit

The following NEW packages will be INSTALLED:

  megahit            bioconda/linux-64::megahit-1.2.9-h8b12597_0

Proceed ([y]/n)?

Where I click yes and the packages is installed. However, after installation I am not able to use the command megahit as this results in the error: megahit: command not found.
I use Ubuntu 18.04.3 LTS and have installed miniconda3.
I reproduced the installation of megahit on a windows computer with a virtual machine using the same Ubuntu version. On my virtual Ubuntu I CAN use the megahit command.
Could anyone please help me solve this problem? I searched the internet but cannot find a good solution on the issue (also not on the github issue page of this package).

Why not just install Anaconda3?

/home/datura-workstation/miniconda3/bin/megahit ? it should be there

Did you run make install after building?

@ATpoint Ah, gotcha. Thanks

In minconda3/bin there is no executable file named megahit. (I checked the folder on my test virtual machine and there the file is indeed there.) That seems to be the problem. I can find the megahit on my workstation in the directory /home/datura-workstation/miniconda3/pkgs/megahit-1.2.9-h8b12597_0/bin. Why is it there? Can I just copy these files to the miniconda/bin?  All files of megahit are placed in this directory while they should be located in the miniconda3/bin directory.

EDIT: I re-installed megahit and now there the megahit files are present in the miniconda3/bin. However, I still receive the error.

EDIT2: after activating conda activate base, megahit became functional. Thanks all for the help!

type conda info --envs.  Look at your list of environments.  Say it's base that has megahit in it, do conda activate base, then just do megahit.

even in the base environment, the command megahit is not found.

I re-installed megahit and now when I activate conda activate base, megahit does work. Thank you!!


text:
How to calculate overlapping genes between two genome annotation versions
I have two annotations of the same genome generated with different annotation pipelines. I want to identify overlapping gene models.
An important feature of this genome is that there are many 'genes within genes', i.e. a genemodel in the intron of another genemodel. Therefore, I only want to count two genemodels as overlapping when their coding sequence exon annotations overlap.
Using something like bedtools intersect it is straightforward to calculate overlap between the gene-level annotations.
However: I am not sure how to select genes as overlapping when only their coding sequence exons (CDS features) overlap.

Why not extract the coordinates of your CDS regions from your bed/gff files, then do the bedtools intersect?

That would still only leave me with the overlapping CDSs features. In the end I want to know the genes. Why not write your comment into an answer?

Short Answer:
In my opinion, my approach would be to pull out the CDS exons and run bedtools on those.
A Few More Details:
When you pull out the exons, make sure that you assign them all IDs if the don't already have them assigned and record which IDs "belong" to which genes. Now when you get exons that overlap, you know that they are coding and you can tie them back to which genes they originate from.

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


text:
How to calculate the number of peaks that are upstream/downstream of some other peaks
I have 3 histone marks,I have used Macs2 for peak calling and diffBind to analyze the peaks. I was wondering if you know any way to calculate the peak numbers of one specific histone mark that occur for example within 2kb of peaks belong to another histone mark.
Much appreciated

You can use windowBed from bedtools with -c if you just want counts. Or use -wa -wb and then use groupBy for any other operations.


text:
merging two/or more bed file into one bed file
I am trying to merge two bed files (more in future) to one. my bed files are something like :
.
I need to merge them in a way to have the shared chromosome location.
Is there a way to do that ?

Convert each of your Excel spreadsheets to tab-delimited text files (how-to).
Install and run Cygwin, if using Windows. If you are using OS X, open up the Terminal app.
Install the BEDOPS toolkit to get sorting and set operation tools (if you're not averse to using them).
Strip out the Microsoft line endings, strip the header, and sort each of the N files:
$ tr -d '\r' < fileN.tsv | tail -n+2 - | sort-bed - > fileN.bed

Merge the files:
$ bedops --merge file1.bed file2.bed ... fileN.bed > answer.bed

The file answer.bed will contain all merged regions from files file1.bed through fileN.bed. Merging will be done on intervals for each unique chromosome name.

Thanks that worked

In my opinion, the easiest way to merge bed files is to use bedtools merge.

To merge multiple BED files into one, you need to do the following:

concatenate them all (bedtools merge only accepts one input BED)
sort the concatenated BED
merge the sorted BED with bedtools merge

Here is a simple pipeline to do this
cat a.bed b.bed | sort -k 1,1 -k2,2n | bedtools merge > merged.bed


text:
How to make Venn diagram for the Macs peak calling output of two data sets?
I have two outputs of macs2 peak calling for two of my data sets. I wanted to plot the Venn diagram to see how many peaks are shared. I mainly work on Unix/Linux. Do you know any way that I can have those diagrams plotted?
Much Appreciated

Venn diagrams can be helpful, but for peak calling they can be somewhat misleading.
This topic is addressed in the DiffBind documentation, in the "Comparison of occupancy and affinity based analyses" section.
Suffice to say, asking "which peaks are shared?" is not the same as asking "is there evidence for differential peaks, and where are they?".
For this reason, performing a differential binding analysis might be a better option.
If you're just looking for Venn diagram plotting software, though, the VennDiagram (note, I have personally contributed to this tool in the past) and UpSet R packages are useful tools.

Check out the tool Intervene. It's exactly designed for this. Github link.

Dear @Cat if you describe Intervene for those of us who are unfamiliar with it?


text:
differential analysis of chip-seq data
I have several sets of chip-seq data. I called the peaks using Macs2. I am pretty new to the field and I will appreciate any help. I wanted to annotate the peaks and see which peaks are shared between them and which ones are different. What is the best way to do this? which package should I use?
Many thanks

For sharing, you can test the proportion of overlapping peaks using bedtools intersect or the find.overlap function in GenomicRanges packages in R.
What makes more sense is perhaps to merge all the peaks together using bedtools merge, and then count the number of reads/fragments in each peak for each sample using featureCounts. The output of this is a matrix, which is very similar to gene expression results. Then you can do some differential tests on this data.

In a similar question, I recommended that directly comparing called peaks can be somewhat misleading. This topic is addressed in the documentation of the DiffBind R package on Bioconductor, in the "Comparison of occupancy and affinity based analyses" section.
But if you're just looking for genome annotation software, the annotatr package on Bioconductor is pretty useful

It depends on how many samples you have.
If you only compare one sample to the second, you can pool and call peaks and then directly calculate RPKM ratios for each peak.
It you have more samples, you may want to use IDR to score differential peaks.
You may also use packages such as Fseq or HOMER to treat one group as ChIP experiment while the other group as control to call peaks.
Alternatively, you could pool all the samples to call peaks, then calculate read counts in each peak for each sample, and then use DESeq2 or other differential tools to calculate significantly different peaks’ significance and effect size.


text:
annotation using ChIPseeker package error
I have differential binding sites object obtain from diffBind (dba.report). I am using the ChIP Seeker package to annotate the peaks but keep getting the following error: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’
I used diffloop package to add chr to the chromosome column but it did not work. Result in the same error.
Please let me know if you have any idea how to fix it. Highly appreciated
Here is my peakset
BD

  seqnames              ranges strand |      Conc Conc_condition Conc_Control      Fold      p-value         FDR
               <Rle>           <IRanges>  <Rle> | <numeric>     <numeric>    <numeric> <numeric>    <numeric>   <numeric>
    [1] NC_000003.12   75668856-75670016      * |  12.27197    12.8243533      0.76175  12.06260 1.43013e-102 1.41852e-97
    [2] NC_000024.10   11320736-11323578      * |  10.58988    11.1419124      0.76175  10.38016 1.15927e-100 5.74930e-96
    [3] NC_000004.12         37431-38634      * |   6.42233     0.0596816      8.06241  -8.00273  6.67673e-98 2.20750e-93
    [4] NC_000017.11   21678164-21686752      * |   9.26507     9.8163416      0.76175   9.05459  3.54723e-95 8.79607e-91
    [5] NC_000022.11   12692637-12694102      * |   9.64945    10.2010195      0.76175   9.43927  2.35838e-87 4.67845e-83
    ...          ...                 ...    ... .       ...           ...          ...       ...          ...         ...
  [226] NC_000004.12 131722954-131744030      * |   11.8156      12.36301     5.337173   7.02584  2.54289e-17 4.78151e-16
  [227] NC_000015.10   68134159-68136364      * |    8.2514       8.80099     0.969455   7.83153  2.28289e-16 3.73532e-15
  [228] NC_000005.10   34189082-34194027      * |   10.4671      11.01445     3.990668   7.02378  1.46745e-15 2.12829e-14
  [229] NC_000018.10       108003-108613      * |   11.2783       5.46414    12.912836  -7.44869  3.68643e-12 3.26794e-11
  [230] NC_000004.12   49709306-49710113      * |   12.0662       6.40350    13.698765  -7.29526  1.23922e-08 6.37694e-0

peakAnno <- annotatePeak(BD, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
> peakAnno <- annotatePeak(BD, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
>> preparing features information...         2020-11-12 01:08:37 PM
>> identifying nearest features...       2020-11-12 01:08:37 PM
>> calculating distance from peak to TSS...  2020-11-12 01:08:38 PM
>> assigning genomic annotation...       2020-11-12 01:08:38 PM
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’```

that was solved by changing seq level:
newnames <- paste0(c("1","2","3","4", "5","6","7", "8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","X","Y","MT"))
names(newnames) <- paste0(c("NC_000001.11","NC_000002.12","NC_000003.12","NC_000004.12", "NC_000005.10","NC_000006.12","NC_000007.14", "NC_000008.11","NC_000009.12","NC_000010.11","NC_000011.10","NC_000012.12","NC_000013.11","NC_000014.9","NC_000015.10","NC_000016.10","NC_000017.11","NC_000018.10","NC_000019.10","NC_000020.11","NC_000021.9","NC_000022.11","NC_000023.11","NC_000024.10","NC_012920.1"))

peaks_with_new_seq_level <- renameSeqlevels(peaks_with_old_seq_level,newnames)```


text:
How to identify specific motif for chip-seq data
For chip-seq data, I had analyzed the samples until motif discovery using both cluster profiler and Homer for each histone mark. I was wondering what is the efficient method to determine if each GO term or identified motifs has a biological meaning or not or if there is any way to filter out the genes/motifs for only specific cell types such as PBMCs, how to prioritize genes? I know random matches might even surpass the biologically relevant matches but right now I am not sure how to narrow down the motifs/genes? Your help would be greatly appreciated.

To add more detail I have extracted PBMCs from infectious patients with Anthrax. I performed chip-seq for 3 histone marks: (H3K27ac, H3K4me3, H3K27me3).I'd like to know if the bacteria is affecting the epigenetic signature of PBMCs?
After peak calling and annotation for each marks, I have more than 2000 genes for each mark and more than 500 motifs called by Homer. I am trying to narrow the genes to those which I expect to see in PBMCs cells and those which have biological meaning. Right now I am struggling to find a way to filter out meaningful genes.


text:
Cannot blast against specific NCBI databases
I am having issues with some prokaryote reference genome databases (exact names : ref_prok_rep_genomes.*), that I downloaded from the NCBI website : https://ftp.ncbi.nlm.nih.gov/blast/db/. Files in this database go from ref_prok_rep_genomes.00 to ref_prok_rep_genomes.12.
In Linux command line, I was able to launch successfully the bash command  blastn -query (my_query) -db ref_prok_re_.genomes.* -out (my_output) with indexes from 00 to 05, but this does not work from indexes 06 to 12 (with the same query and output files). The exact same thing occurs with blastdbcmd -info -db ref_prok_rep_genomes.*.
Each time the error displayed is the same:
BLAST Database error: Input db vol does not match lmdb vol.
All the ref_prok_rep_genomes files are in the same folder. I downloaded them using wget and extracted them using tar -zxvf. All ref_prok_rep_genomes.* databases (indexes 00 to 12) have nhr, nin, nnd, nni, nog and nsq files. The folder also contains ref_prok_rep_genomes files (without numbers), with  nal, ndb, nos, not, ntf and nto extensions respectively, as well as taxdb.btd and taxdb.bti files.
Do you guys have an idea of what's going on?
Thank you very much for your time :)

Did you checked the integrity of the downloaded files? wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/ref_prok_rep_genomes.??.tar.gz.md5 && md5sum -c *md5

I just did. With any of the md5 files, it gave me
ref_prok_rep_genomes.01.tar.gz: OK
ref_prok_rep_genomes.02.tar.gz: OK
ref_prok_rep_genomes.03.tar.gz: OK
ref_prok_rep_genomes.04.tar.gz: OK
ref_prok_rep_genomes.05.tar.gz: OK
ref_prok_rep_genomes.06.tar.gz: OK
ref_prok_rep_genomes.07.tar.gz: OK
ref_prok_rep_genomes.08.tar.gz: OK
ref_prok_rep_genomes.09.tar.gz: OK
ref_prok_rep_genomes.10.tar.gz: OK
ref_prok_rep_genomes.11.tar.gz: OK
ref_prok_rep_genomes.12.tar.gz: OK```

Actually I understood what happened. NCBI's new version 5 databases use virtual LMDB memory. They give issues with the previous versions of ncbi-blast+.
I managed to fix the problem by (1) updating ncbi-blast+ (currently it is the 2.10.0 version), (2) deleting my previous databases, (3) downloading them again using update_blastdb.pl --decompress <database_name> instead of wget


text:
Compressing the human genome to few megabytes
Multiple sources (see for instance this or this) discuss how genetic data will have scalability problems, given the huge file size of the human genome. The most straightforward encoding (see here) of the human genome requires about 700 Mb.
I came across this paper claiming to be able to store the human genome in about 4 Mb, having a reference genome and exploiting the fact that all human genomes are mostly equal. This paper is much older than the the other references discussing the scalability problems. Why is this technique not widely used?

https://bioinformatics.stackexchange.com/questions/12824/theoretical-limit-of-human-genome-compression?noredirect=1&lq=1

Answer from @alex-reynolds, converted from comments:
As with any compression scheme, if you have to store a reference genome, along with the 4 Mb "diff", is it really compressed?
http://mattmahoney.net/dc/dce.html
For sequencing data, you might also look at the CRAM format, which also uses diffs against a reference genome to store reduced datasets. Labs are switching from BAM to CRAM to manage storage of increased numbers of sequencing datasets, with associated changes to pipelines and toolkits to use this format. For performance reasons, BAM may be preferred for analysis, and CRAM for long-term storage.
https://gatk.broadinstitute.org/hc/en-us/articles/360035890791-SAM-or-BAM-or-CRAM-Mapped-sequence-data-formats
https://www.ga4gh.org/news/cram-compression-for-genomics/

I saw this paper a few weeks back that claims to offer 4x more compression than typical formats for both fastqs and BAM/CRAM files.
Genozip - A Universal Extensible Genomic Data Compressor
https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btab102/6135077


text:
What RNA-Seq expression value would be closest to Microarray equivalent?
I know this question may seem strange.
I'm using Spearman correlation between gene expression profiles for various reasons (I won't go into details here). As a result, I often compare RNA-Seq and Microarray samples.
For preliminary analysis, I usually grab what version of the data is easily accessible (RPKM, FPKM,...), but I'd like to dig a bit more.
Intuitively, I'd think a value similar to RPKM or such would make more sense than raw counts, which is why I usually convert raw counts to RPKM (I know they're considered obsolete for statistical analysis of the data, what I'm doing relies on the correlation)
So I ask what RNA-Seq gene expression value would be closest to a microarray equivalent ?
i.e. What would (theoretically) maximize correlation between gene expression of the  same sample profiled with RNA-seq and microarray. RPKM ? FPKM ? TPM ? etc.

“I know they're considered obsolete for statistical analysis of the data” — No, they’re considered obsolete for every purpose. Just don’t use them. Read https://bioinformatics.stackexchange.com/a/69/29.

I did a comparison of cDNA count data against microarray data that was published a few years ago:

For comparisons to published data (Fig. S2; Miller et al., 2012), a generalized linear model was fitted to the relationship between log-transformed microarray and VSTPk expression levels obtained from the ImmGen Project database, and was used to transform the microarray data into values comparable to our VSTPks.

I found that the Variance-Stabilizing Transformation that was carried out by DESeq2 was close to what I wanted, but there still seemed to be a length-based bias to the reads. I corrected this by dividing by the length of the longest gene isoform in kilobases (creating something that I called VSTPk).
After doing this, there were range differences between the microarray and cDNA data, so I did an additional linear transformation to get the data fitting as close as possible.

I think it is very hard to say which are the closest because they are not really comparable. But since you are using Spearman correlation, I guess RPKM, FPKM, and TPM do not change the order of gene expression levels. You might also want to normalize RNA-seq and microarray data so that they are more comparable.


text:
What do "PnBS" and "PnMS" mean in FCS flowCore parameters?
I'm using the flowCore packages to analyze cytometry data, and I am trying to understand all the parameters stored in the FCS file.
I've figured out most of the elements, but despite a thorough search I haven't been able to find documentation on some of the parameters.
For each channel, there are 9 associated parameters. Here's P1 and its associated values from one of my files :

$P1N    FSC-A
$P1R    262144
$P1B    32
$P1E    0,0
$P1V    225
$P1G    1
$P1DISPLAY  LIN
$P1BS   0
$P1MS   0

I've figured out the following from various sources (notably the FCS data standards, given here) or from the fields themselves:

$PnN Short name for parameter n.
$PnR Range for parameter number n.
$PnB Number of bits reserved for parameter number n.
$PnE Amplification type for parameter n.
$PnV Detector voltage for parameter n.
$PnG Amplifier gain used for acquisition of parameter n.
$PnDISPLAY  the transformation to apply when displaying the data (e.g. LIN, LOG)

However, these 2 don't appear in the file standard description, and flowCore does not have documentation on the parameters either.

$PnBS
$PnMS

I've seen $PnBS set at -1, 0, 48, 57, 54, 52 or 25 for various channels in my data, while $PnMS always seems to be at 0.
Does someone know what it is, or where to look for a definition?
If that helps, my FCS files were generated through BD FACSDiva Software Version 9.0.1.

As this is a package from Bioconductor I would recommend to ask on their site: support.bioconductor.org. Probably you'll get a helpful answer, but it would be nice to know where do you find these parameters

I found my answer in an old mailing list. Credit to Josef Spidlen who authored the mail.
Short answer: both parameters are proprietary keywords from the BD FACSDiva software for storing scaling values.
Long answer from the mailing list:

PnBS Indicates the automatically set Biexponential scale value for Instrument  Parameter (-1 indicates an invalid scale value). When
importing an FCS  file, if the FCS file has valid Biexponential
scales, they should be  applied to the data. If the FCS file doesn't
have valid biexponential  scales then the scales are calculated from
the entire data. If  biexponential scales exist, then if AUTOBS is
set, the PnBS values are  used and applied to the data. If AUTOBS is
not set, the PnMS values are  used and applied to the data. The
information represented by this keyword becomes part of a Tube's
attributes, and does not appear in the Tube's list of keywords. It can
be viewed in the Biexponential inspector' when the auto scaling option
is selected.
PnMS Indicates the manually set Biexponential scale value for Instrument  Parameter (-1 indicates an invalid scale value). When
importing an FCS  file, if the FCS file has valid Biexponential
scales, they should be  applied to the data. If the FCS file doesn't
have valid biexponential  scales then the scales are calculated from
the entire data. If  biexponential scales exist, then if AUTOBS is
set, the PnBS values are  used and applied to the data. If AUTOBS is
not set, the PnMS values are  used and applied to the data. The
information represented by this keyword becomes part of a Tube's
attributes, and does not appear in the Tube's list of keywords. It can
be viewed in the Biexponential inspector' when the manual scaling
option  is selected.


text:
How can I extract all known mutations of my BAM (or SNP/INDL files)?
I am using a Genome Explorer tool to see all the mutations on my own DNA.
My particular interest is on listing the variants and get their names/ids.
Here are a few screenshots:

You can see that they tell you the genomic position, the allele frequency of the mutation, the name of the disease, the Category, etc etc...
However, the filters on this tool work very badly so I cannot filter all my mutations only...
So I was trying to find out an alternative tool to get this data in such a nice way this tool provides. I tried IGV but they just provide raw readings, I didnt find a way to get a list of variants the same way this Genome Explorer does.
The closest thing I've found was Geneious:

but it crashes on my computer and is too complex to use (I wasn't able to reach that screen above, from my BAM file). And it also eats a lot of RAM that I don't have.
Are there any alternatives? Even paying?

If you have a BAM file, you could easily run a variant caller analysis software such as freebayes, though you would need a good amount of computing power like a cluster or a dedicated node with good memory. Freebayes outputs a vcf file that has all the variants listed. Then you can load it into the UCSC or Ensembl human genome browsers and look at what genes they are.

Hi! Thank you for sharing! I already have these files as well actually (.snp.vcf.gz, .indel.vcf.gz, .cnv.vcf.gz and their .tbi's). Can you detail me a bit more on the part where I can load these into some tools to list all mutations known to mankind at current 2020 (and their related disease/classification)? I am exploring both websites you just mentioned but feeling a bit lost. Thanks!

I think this one does what I want (http://www.ensembl.org/Homo_sapiens/Tools/VEP), but is giving me "bad gateway".
Any chance this is an open source project?

Found it: https://github.com/Ensembl/ensembl-vep

Yes vcf, a discussion on vcf algorithms would be cool


text:
Use a specific python.exe instance with PyMOL/Windows?
I want to use RDKit with a python script in PyMOL. RDKit requires a conda environment, so I think I need to be able to specify that the instance of python that I want to use would be in that environment, instead of the default for PyMOL, so that I can run it through PyMOL. The script I made is working for me through PyCharm when I test it, by using the RDkit environment "rdkit-env" as my python interpreter, with all of the additional modules I need installed there as well.
I would like to run it through PyMOL so I can add some PyMOL specific functionality as well. I've tried using shebang to specify the RDkit environment that I made, but that does not seem to be working for me with Windows.
Could someone please help me figure out what other approach might work for this? I'm new to this area, so I would appreciate any insights. Thank you!

So it is really hard to make a given python instance run a library from a different instance. It is not just one dist-packages library folder, but multiple that you would need to link, and all that assuming they use the same version of python.
So the best bet is to install a new instance of pymol the python module.
The module can be installed in a python that uses conda with
conda install -c schrodinger pymol

This is the incentive version of PyMOL 2 and does not require a licence.
You installed rdkit, so I am assuming you used conda as conda install -c rdkit rdkit. But just in case you went through the pain of compiling rdkit for a given python instance without conda and want the same for PyMOL. There are guides of instructions required to compile open source. Such as https://pymolwiki.org/index.php/Windows_Install. On a barebone OS like a Pi's Raspsian, I had to do a lot more than claimed say sudo apt-get install python3-dev libglm-dev freeglut3-dev libglew-dev libpng12-dev libfreetype6-dev build-essential libxml++3.6-dev and link a few things manually. So just like for rdkit, conda make life a breeze.
PS. RDKit and PyMOL? Add PyRosetta and you just my favourite modules!


text:
Is there a way to save a full pdb file in PYMOL after uploading biological unit?
pdb code: 2obt
The pdb file has a coordinates of a monomer.
I have used following commands to generate an assembly, or in other words to upload one more biological unit so that the protein would be a dimer:
set assembly, ""
fetch 2obt, asu, async=0

set assembly, "1"
fetch 2obt, assembly1, async=0

After generating the functional dimer in  pymol under name assembly1, I have an important issue that I cannot save a dimer structure in pdb file with dimer coordinates. Does anybody has an idea how to go around it or is it even possible to do that?

can we says that "after uploading the biological unit is not the right wording" ?? --> When loading mmCIF structure files, the assembly setting controls whether PyMOL loads the asymmetric unit (assembly="") or a biological unit (e.g. assembly="1"). Assembly identifiers in mmCIF files can be arbitrary strings, but usually are numeric and most files define at least assembly "1".  I guess the biological unit is created applying same symmetry operator on the asymmetric unit (or am I wrong ??). Ok PymolWiki says loading, actually loads .....

So on the GUI form of PyMOL on the right hand side you have a list of what is ambiguously called objects in PyMOL (models elsewhere), which can be enabled or disabled by clicking on their name.
The name also has a 1/1 in green or similar. This indicates the state. In your case you have 2 states.
To split a multistate object use
split_state prot_A

To combine two objects one uses the command
create new_name, prot_A or prot_B

Now, this has two caveats which if violated result in a nasty piece of modern art:

if the chains, segment ids and residues are the same
if you had two overlapping macromolecules

To deal with the first, change the chain id on one of them.
alter prot_B, chain='B'
   sort
Segment identifier is a weird thing that most comp-biochem tools dislike, so I like to remove that
alter prot_B, segi=''
   sort
Now when you go to save your new_name object, in the popup of the GUI you get a dropdown menu defaulted to enabled (visible objects that have not be clicked to be disabled). There you can choose new_name.
Aligning two
If someone that is in not the OP author has a template protein and needs to make a dimer of another based on that, the process is the same.
create twin, prot_B
align twin, prot_A and chain B
alter twin, chain='B'
sort

SDF files
If someone else has a multicompound sdf where the state is the header name, I strongly advise not to use split_state. Yes, the python pymol module (not the app) handles states in a really weird way, but states can be selected with state 1 selection algebra and state is an argument for most commands that do not accept state selector.


text:
Calculating polarity and hydrophobicity of a ligand/small molecule
I have lot of SMILES IDs of the ligands. What I am interested in is the polarity and hydrophobicity score of each of the ligands.
In other words, does ligand have more polar groups and can it make more polar interactions or is it more hydrophobic?
Does anybody know the way/software to look at it?
Best

You can calculate their polar surface area or logo using rdkit or openbabel. Other molecular properties may be available as well.

I see, you mean like TPSA of the molecule? Also, is there some cutoff that can be used there, in a way like TPSA of 100-150+ is considered that molecule is polar, while TPSA<100 is considered as hydrophobic molecule?

I am not sure what are the threshold or if there is a threshold or you look at literature


text:
Visualizing protein-ligand interactions with LIGPLOT or any other suitable software?
I have used DOCK software to dock 10 ligands into the protein structure.
I want to generate a pdb file that contains all 10 binding modes of the ligands together with the protein. So all together in one pdb file.
Now, what I want to know is when I give this file to LIGPLOT, is it possible to generate something like:
%polar interactions XX
%non-polar interactions XX
%hydrogen interactions XX
Is there a way to do it with LIGPLOT, if now is there some other software that can be recommended for this type of assessment?
I have been having hard time, trying to find a way or software that can perform similar type of calculation and give me similar type of output that I want.
Best

"binding modes" --> "poses".
You will need to do it programmatically. This post for example use Python pymol for querying atoms. Can you use Python, Perl or Julia?

Yeah, I assumed so, cus it tends to be a bit more complicated. Unfortunately, I cant say for my self that I know any of the three languages :(

You can visualize with Discovery Studio Visualizer and Schrodinger.
Arpeggio and PLIP can generate output files of the different interactions.

Hey, thanks! already found that PLIP can do it and it works pretty well so I went for it.


text:
How can I convert codon coordinates to genomic position?
I am looking for a given mutation in IGV, which accepts coordinates in the form of
chr<X>:<Y>

Where X is the chromosome number and Y is the base location in that chromosome.
Eg. chr6:26091179 which has the following attributes:
Gene: HFE|LOC108783645
Marker: rs1799945
Position: chr6:26091179

Additionally, an interval can be written like
chr<X>:<Y>-<Z>

So basically from base located at Y to Z.
My questions are:
1) How do I convert HFE|LOC108783645 to chr6:26091179 (and vice-versa)?
2) How do I convert rs1799945 to chr6:26091179 (and vice-versa)?
3) On a different example, how do I convert something like codon 109 of exon 3 of beta globin gene (GTG->TTG) to chr<X>:<Y> (or chr<X>:<Y>-<Z>)?
Why are there so many ways to specify positions in genes? For me the most clear way is chr<X>:<Y>, all the others leave me guessing...

You should look at HGVS conventions used to refer to positions with respect to reference sequences. Remember, not even chr:pos notations are independent of reference sequence versions. chr6:26091179 in GRCh38 is not the same in GRCh37.

HFE|LOC108783645 is https://www.ncbi.nlm.nih.gov/gene/108783645. You can use biomaRt to go from the LOCxxxx ID to location. Searching by HFE will give you a different location though, so you might have to account for such eventualities. It is a long and often manually intensive process.
You should be able to use entrez-direct/eutils to go from rsID to position
Unless there's some sort of natural language processing API, codon 109 of exon 3 of beta globin gene cannot be parsed by existing APIs (Unless it's always codon X of exon Y of Z, in which case you could use grep too get the components. Even so, going from exon + codon + codon base change to exact base location will require a lot of wrangling.


text:
Problem with classification model of genomic data: every machine learing model predicts wrongly almost always the same subset of dataset
First of all, I'd like to apologize for any spelling or grammar mistakes.
I'm having a problem using R for a classification problem. My dataset contains ~300.000 genomic data, and the features are DNA-related features (number of dinucleotides, number of trinucleotides, the CG Content, and some more). In conclusion, I have a dataset of 300.000 rows and 84 columns (columns= features). The 84th feature is basically the classification variable (there are two classes: class 1 and class 2).
I split the dataset in train set and test set and fit the train set to create a model using different methods, like Random Forest, XGBoost, Adaboost, etc. and check the accuracy using the test set.
In every model, the accuracy of the prediction on the test set is always around 86-87%. What I noticed is that, in the test set, for EVERY MODEL, almost the same number of rows (data) is predicted right (the model predict its true class), and always almost the same data of rows (data) is predicted wrong (the model predict wrong class. For example its true class is 1 but the model predicts 2).
If I take the set of data that are predicted right for every model I've created (250.000 data from the 300.000), and build a new model out of them (splitting them again in a new train and test set), I have an accuracy of 98%.
That means that the other 50.000 data (from the 300.000), add an error at every model I'm trying to create.
But I can't just ignore them. I thought that it would be a way to create two differents models, but doesn't know how to figure this out.
Then I thought that there could be a way to find which variables from the two subsets have strong differences, but even if I figure this out, I don't know what's the way for example to scale the original dataset to prevent those differences from keeping the accuracy in 87% in every model (using the original dataset of the 300.000 rows).
In conclusion, I have a dataset that pretends like "two datasets". It consists of two subsets: the first one is always predicted right from any model (machine learning algorithm) I'm trying, and the second one is always predicted wrong. How could I build a strong model taking into account the differences?
I also took a look at the box plots of the features. The most important feature for every model is called "CG" (number of dinucleotide of cytosines and guanines in DNA). In the following image you can see the boxplot of CG feature separating the original dataset in two datasets (T = for the dataset that is predicted right in every model) and (F = for the dataset that is predicted wrong in every model). 1 stands for data of class 1 and 2 stands for data of class 2. So, T1 means data from the right predicted dataset of class 1, T2 means data from the right predicted dataset of classes 2, T means data from the right predicted dataset for both class 1 and 2, etc.

As I'm a new beginner in R and in machine learning, I hope you do understand my problem. Thanks sincerely.

It needs some thought. Whats your training/test split, 80:20? Any parameter splits?

The persistent 87-88% accuracy between models is a bit weird, but its an interesting problem

I always used a train set with same amount of class 1 and class 2. for example, a train set of 50.000 data (25.000 of class 1 and 25.000 of class 2).

Okay thats the classificaiton (I don't know what "Class I" and "Class 2" represent, but I get the idea in ML, whatever "Class" is, is the thing you are looking at. The question was the split with testing, or parameterisation. The reliability of the accuracy score across very different models suggests something is not right, but I honestly don't know what.

Are your classes balanced, meaning that do you have roughly the same number of the two classes? Looks like your models are learning a particular set of samples very well at the expense of learning others. Cross-validation could help here.

Not sure it is related but do you shuffle your datasets before training ? I had weird patterns of valid accuracy from epoch to epoch when training a neural network once, it was because I forgot to shuffle my datasets.

You could be seeing the effect of SNPs on haplotype-driven conditions (i.e. where a combination of multiple variants predicts the class, rather than a single variant).
Do you have any method to determine the age of particular variants (e.g. via an LD plot to find inconsistencies between variants and the surrounding block structure), so that you can see if that correlates with the likelihood of error?

No, unforetunately no.


text:
Weighted Network Analysis Book by Steve Horvath, taking the power of a matrix
In Horvath's Weighted network analysis, the author mentions that M^b,  where b is a real number and M is a matrix, denotes the element-wise power. This seems odd to me. The way I was taught, M^b would mean multiplying the matrix M, b times. Should element-wise power not be M^(ob) where where o denotes the Hadamard power?

You CAN calculate M^b for b being float number, complex number or even another matrix.
As well as you can calculate exp(M), log(M), sin(M) whatever... That is very standard for math people.
Let me give two ways how you can think of it.
1) M^b = exp(b log(M)), so you need to define what is exp and log for matrices.
That can be done using power series: exp(M) = 1 + M + M^2/2! + M^3/3! + ...
log(1+M) = M - M^2/2 + M^3/3 - ...
2) More general way to calculate any function of matrix f(M) , looks like that:
if your matrix is diagonal, then it is true that: f(diag(d_i)) = diag( f(d_i)  ) .
So the question how to reduce any matrix to diagonal one.
That can be done for almost all matrices by conjugation M = C D C^{-1}.
Then f(M) = C f(D) C^{-1}.


text:
Tools for creating gene coexpression networks
What are some good tools (in Python if possible) for creating gene co-expression networks from RNA seq data?

it depends on the size of your data. If you have enough RAM available, then you may directly compute a correlation-based distance matrix using the cor() function (with method="spearman" or method="pearson"). But you may also give a try to Mutual Information distance using the Parmigene (https://cran.r-project.org/web/packages/parmigene/index.html) R package or partial correlations with the Corpcor (https://cran.r-project.org/web/packages/corpcor/index.html) R package.
If you have a large expression matrix, remember to use the excellent data.table package and its fread() function.

After looking around, it seems like the Camoco library is suitable.

It would be nice to hear from the R crowd


text:
xtract from eutilities not exctracting the information I need
I'm running the following command which uses the NCBI e-utilities:
esearch -db bioproject -query "PRJNA198476 [PRJA]" | esummary
This returns:
<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE DocumentSummarySet PUBLIC "-//NLM//DTD esummary bioproject 20140903//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20140903/esummary_bioproject.dtd">

<DocumentSummarySet status="OK">
<DbBuild>Build211111-1430.1</DbBuild>

<DocumentSummary><Id>198476</Id>
    <TaxId>6239</TaxId>
    <Project_Id>198476</Project_Id>
    <Project_Acc>PRJNA198476</Project_Acc>
    <Project_Type>Primary submission</Project_Type>
    <Project_Data_Type>Transcriptome or Gene expression</Project_Data_Type>
    <Sort_By_ProjectType>476331</Sort_By_ProjectType>
    <Sort_By_DataType>472932</Sort_By_DataType>
    <Sort_By_Organism>51608</Sort_By_Organism>
    <Project_Subtype></Project_Subtype>
    <Project_Target_Scope>Multiisolate</Project_Target_Scope>
    <Project_Target_Material>Transcriptome</Project_Target_Material>
    <Project_Target_Capture>Whole</Project_Target_Capture>
    <Project_MethodType>Sequencing</Project_MethodType>
    <Project_Method></Project_Method>
    <Project_Objectives_List>
        <Project_Objectives_Struct>
            <Project_ObjectivesType>Expression</Project_ObjectivesType>
            <Project_Objectives></Project_Objectives>
        </Project_Objectives_Struct>
    </Project_Objectives_List>
    <Registration_Date>2013/04/22 00:00</Registration_Date>
    <Project_Name>Caenorhabditis elegans</Project_Name>
    <Project_Title>Deep sequencing of endogenous mRNA from Caenorhabditis elegans in the presence and absence of arsenite</Project_Title>
    <Project_Description>Background: Arsenite is one of the most toxic chemical substances known and is assumed to exert detrimental effects on viability even at lowest concentrations. By contrast and unlike higher concentrations, we here find that exposure to low-dose arsenite promotes growth of cultured mammalian cells. In the nematode C. elegans, low-dose arsenite promotes resistance against thermal and chemical stressors, and extends lifespan of this metazoan, whereas higher concentrations reduce longevity. While arsenite causes a transient increase in reactive oxygen species (ROS) levels in C. elegans, co-exposure to ROS scavengers prevents the lifespan-extending capabilities of arsenite, indicating that transiently increased ROS levels act as transducers of arsenite effects on lifespan, a process known as mitohormesis. The RNA-seq data comprises 2 biological replicates for worms exposed to 100nM Arsenite 48h after L4 and 2 biological replicates of the same age as controls Jena Centre for Systems Biology of Ageing - JenAge (www.jenage.de) Overall design: 4 samples: 2 mRNA profiles of C.elegans 48h after L4 exposed to Arsenite; 2 mRNA profiles of C.elegans 48h after L4 as controls (H20). The N2 wild type (var. Bristol) strain was used.</Project_Description>
    <Keyword></Keyword>
    <Relevance_Agricultural></Relevance_Agricultural>
    <Relevance_Medical></Relevance_Medical>
    <Relevance_Industrial></Relevance_Industrial>
    <Relevance_Environmental></Relevance_Environmental>
    <Relevance_Evolution></Relevance_Evolution>
    <Relevance_Model>Yes</Relevance_Model>
    <Relevance_Other></Relevance_Other>
    <Organism_Name>Caenorhabditis elegans</Organism_Name>
    <Organism_Strain></Organism_Strain>
    <Organism_Label></Organism_Label>
    <Sequencing_Status>SRA/Trace</Sequencing_Status>
    <Submitter_Organization>Leibniz Institute for Age Research - Fritz Lipmann Institute</Submitter_Organization>
    <Submitter_Organization_List>
        <string>Leibniz Institute for Age Research - Fritz Lipmann Institute</string>
    </Submitter_Organization_List>
    <Supergroup>Eukaryotes</Supergroup>
</DocumentSummary>

</DocumentSummarySet>

I'm trying to extract the Submitter_Organization using:
esearch -db bioproject -query "PRJNA198476 [PRJA]" | esummary | xtract -pattern DocummentSummary -element Submitter_Organization
However, this returns nothing, and I can't figure out why

You have a typo in DocumentSummary.

@vkkodali thank you that was it


text:
How can I programmatically retrieve all alternative gene names from Wormbase?
I'm trying to get all C. elegans genes from Wormbase programmatically, as well as all of their alternative gene names. To do this, I'm using intermine.
I was able to obtain all genes, but not ALL alternative names. The code below retrieves the genes' primary identifiers, and their "secondary identifiers."
from intermine.webservice import Service
service = Service("http://intermine.wormbase.org/tools/wormmine/service")

query=service.new_query("Gene")
# 6239 is the taxonId for C. elegans
query.add_constraint("organism.taxonId","=","6239")
query.select("*")

all_wormbase_genes_symbol = list()
all_wormbase_genes_sID = list()

for row in query.rows(start=0):
    # get primary symbol/ID
    all_wormbase_genes_symbol.append(row["symbol"])
    # get secondary ID
    all_wormbase_genes_sID.append(row["secondaryIdentifier"])

I tried to look for other keys in row that might contain this information by first retrieving the list of all possible keys:
k = set()
for row in query.rows(start=0):
    k.update(row.keys())

but none have what I want:
{'Gene.automatedDescription',
 'Gene.biotype',
 'Gene.briefDescription',
 'Gene.description',
 'Gene.geneticPosition',
 'Gene.id',
 'Gene.lastUpdated',
 'Gene.length',
 'Gene.name',
 'Gene.operon',
 'Gene.primaryIdentifier',
 'Gene.score',
 'Gene.scoreType',
 'Gene.secondaryIdentifier',
 'Gene.symbol'}

I don't know via intermine, but you can download a file with this information from the FTP at /pub/wormbase/releases/WS287/species/c_elegans/PRJNA13758/annotation/c_elegans.PRJNA13758.WS287.geneOtherIDs.txt.gz (changing to the appropriate WSxxx version)


text:
Speeding up network creation in WGCNA
I am following the tutorials on the WGCNA website for creating a coexpression network. It is taking way too long. Is there a way to speed up execution? All I need is to create the co-expression graph. I will then export to Python and do my analysis there.
In particular, the step "calculating module eigengenes" is taking a lot of time. I do not really care about the modules. All I want is the graph with the genes and the weighted edges.

You can use the "manual network construction" approach. Calculate TOM, then use exportNetworkToCytoscape or exportNetworkToVisANT to export the network (or simply save the TOM matrix as a table and load it in Python).

I tried doing exportNetworkToCytoscape(TOM), but I get Error in exportNetworkToCytoscape(TOM) :    Cannot determine node names: nodeNames is NULL and adjMat has no dimnames. How can you let it just use the names that were used when creating the graph?

Supply the names in the argument nodeNames.

So I did as suggested, I ran the command exportNetworkToCytoscape(TOM, nodeNames = names(datExpr), edgeFile = "edges.csv", nodeFile = "nodes.csv") but I'm curious about the result. edges.csv seems to have been populated correctly (though I'm not 100% sure, there are 5237 edges, which sounds reasonable), but nodes.csv not so much. In particular, the nodes.csv file has only 151 node entries, but I'm not sure why, since in the tutorial there are 3600 genes and so I expected 3600 entries.

If you want the entire network, you also need to set argument threshold to zero. otherwise you will only get edges whose weight are above the (default) value of threshold which is 0.5. Just be prepared to see large edge file.


text:
Importing a WGCNA co-expression network into a networkX graph in Python
Is there a way to generate a WGCNA co-expression graph in R, and then to import that graph into Python such that I can represent it using a networkX object?

You can use exportNetworkToCytoscape or exportNetworkToVisANT to generate and egde list file with the two nodes at the ends and weight for each edge. Hopefully there's a way to import that into networkX.


text:
What is the best way to get a large number of RNA seq data from SRA in Python without being denied access
I have the following code to download data from the SRA using multithreading in Python. After running this a few times now (for testing purposes), I keep getting denied access to the data. Not sure how to fix this. In particular, cating the output files gives:
<?xml version="1.0" encoding="UTF-8"?><Error>AccessDenied<Message>Access Denied</Message><RequestId>cb0f0e98-cafb-1dd7-9b7b-d8c49756ec52</RequestId><HostId>QeDGVwBXYp61J0B4_OUTn7UsEsiQEec0n18DAeR0kaE</HostId></Error>

Here is my code:
import threading
import requests
import time

start = time.perf_counter()

class MyThread(threading.Thread):
    def __init__(self, url):
        threading.Thread.__init__(self)
        self.url = url
        self.result = None
        self.filename = url.split('/')[-1]
    def run(self):
        res = requests.get(self.url)
        with open(self.filename, 'wb') as f:
            f.write(res.content)
urls = [
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.2',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.2']

threads = [MyThread(url, ) for url in urls]
for thread in threads:
    thread.start()
for thread in threads:
    thread.join()

finish = time.perf_counter()

print(f'Finished in {round(finish-start, 2)} second(s)')

I once needed to download more than 1000 files from SRA, and I found a blog post that is very useful (https://reneshbedre.github.io/blog/fqutil.html). The idea is to use aspera to fetch SRA files. This is fast becasue of aspera. Then use fasterq-dump to convert SRA file to fastq. fasterq-dump is faster than fastq-dump partly because it does not compress the output file. I think compressing is the rate-limiting step here. After this, you can use pigz, which is a multi-threading version of gzip to compress the fastq files. Note that if you want to download files to a computing cluster, you can submit jobs for fasterq-dump and pigz. This makes things much faster. You cannot submit jobs for fetching SRA, because you need to connect to internet.

Thanks for sharing this. So I tried running prefetch with aspera as explained in the blog, but I get unrecognized option when I include the '-a' option. I do have aspera installed. Did you have this issue by any chance?

seems like aspera is no longer supported https://github.com/ncbi/sra-tools/issues/255

however, prefetch is still very fast, thanks for sharing this, update your answer to reflect this and I'll accept it!


text:
Programmatically retrieve Accession List from NCBI
I would like to programmatically retrieve the file "Accession List" from an NCBI page.
I'm running firefox, so I'm clicking the file and clicking Ctrl + Shift + E to try and figure out what the file endpoint is. But I can't seem to figure it out (I haven't done this much, I just remember doing it once for a previous project).
When I do this, I see two links popping up on my screen, as shown below:

I see that the first one has a field (I'm talking about the panel on the right) called "acc", which varies from one data dump to the next, so I'm presuming this is the link I need to use to retrieve the file (though I'm not sure), and just change the "acc" field depending on the data I want. Assuming this is right, I'm not sure how to go from here.
I've included the text on the right below for clarity.
GET
    https://www.google-analytics.com/collect?v=1&_v=j83&a=1402353855&t=event&_s=4&dl=https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA491191&o=acc_s%3Aa&ul=en-us&de=UTF-8&dt=Run Selector :: NCBI&sd=24-bit&sr=1280x720&vp=1232x601&je=0&ec=Link Click&ea=Metadata  &el=NONE&_u=SACAAEAB~&jid=&gjid=&cid=1350895533.1590937485&tid=UA-68648221-3&_gid=391407234.1594002343&cd1=run_selector&cd2=sra&cd3=&cd5=NO&cd6=run_selector:sra:&cd8=&cd9=&cd10=&cd11=&cd12=&cd13=&cd14=&cd15=&cd16=&cg1=run_selector&cg2=sra&cg3=&cg4=run_selector:sra:&z=508671370

Which accessions are you trying to get to? In the page you refer to, there is an option to download the list of accessions - one SRR run accession per line - by clicking the 'Accessions' button or the entire table by clicking the 'Metadata' button as highlighted in the screenshot below:

You can get the same data from the command line using the Entrez Direct. The command for that would be something like:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo

to download the table in CSV format.


text:
Understanding the "relativeness" in RNA seq experiments
Let's say we have two samples, sample 1 and sample 2, and we're running an RNA seq experiment. Consider genes A and B for both. For simplicity, let's say that when we run the RNA seq experiment, we get the following results:
gene A sample 1 count: 2
gene B sample 1 count: 4
gene A sample 2 count: 2000
gene B sample 2 count: 4000
Based on my understanding, RNA seq gives relative results. So I guess in this guess, how can you tell whether either one of the following are true:
1- You sequenced sample 2 a thousand fold more compared to sample 1
2- Sample 2 really has a thousand fold more expression in both genes A and B

We cannot be sure about this. For all we know, they could have been sequenced to the same depth and there is a biological difference that drives the 1000x (ridiculously large) fold change.
See above

The only thing you can say for sure is that gene B is expressed twice as much as gene A in samples 1 and 2. That's what "relative" means - relative expression of genes within a sample.

thanks this was helpful. just want to see what others have to say as well before I accept the answer

Oh of course. I’m no expert at this, but this is what I’ve encountered in our experiments.

The process of RNA-seq is a sampling process of all the RNA molecules. We cannot make reliable comparisons before we know how deep we sample. And since we know it is a sample from a population, the absolute value does not make sense.

So can we know how deep we sample? And knowing that, what would that allow us to do exactly? And how does this tie to sequencing "depth"?

This is usually the sequencing depth when I think about this problem.


text:
How to tell if library generation for RNA-seq experiments are stranded vs not stranded
I would like to use RNA-seq data from the NIH GEO database. How can I tell if library generation for the RNA-seq experiments are stranded vs not stranded?

If you already have a reference transcriptome assembly, then you may run the Salmon tool for quantification using your RNA-seq samples. It has an automated library type inference function (https://salmon.readthedocs.io/en/latest/salmon.html#what-s-this-libtype). It might help you to determine whether the libs are stranded or not.

Unfortunately, the stranded nature of cDNA sequencing is not typically included in the metadata attached to a project, which means that downloading a small subset of the data and checking for strand consistency (as suggested by the other answer) is often the quickest way to work that out (when that information isn't included in the methods section of the associated paper). This knowledge is not particularly useful for casually browsing projects to find a stranded protocol that fits a particular study specification.
Sometimes you get lucky, but not often. GEO RNASeq datasets usually have links to the SRA accessions, and those accessions usually include annotation about the library type. For example, consider this dataset on GEO:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE151535
It links to SRP265460 on SRA. Following the link for the first run for this project, I get to here:
https://www.ncbi.nlm.nih.gov/sra/SRX8436504[accn]
This includes the following metadata about the run:

Instrument: NextSeq 500
Strategy: RNA-Seq
Source: TRANSCRIPTOMIC
Selection: cDNA
Layout: SINGLE
Construction protocol: Replicates of oropharyngeal explants at different stages of development (37/38, 39, 41) were immersed in RNALater. Samples were homogenized with a Dounce homogenizer in RLT buffer containing ß-mercaptoethanol (Qiagen). The homogenate was then processed using an RNAEasy kit (Qiagen, Hilden, Germany). Replicate samples for each stage were analyzed for RNA quantity and quality on a Nanodrop 8000 spectrophotometer. RNA libraries were prepared for sequencing using standard Illumina protocols

But, as you can see, there's nothing about whether its stranded or not: Illumina has both stranded and unstranded "standard" protocols.


text:
Files for paired end RNA sequencing
I am looking at the videos at a DIY Transcriptomics course and the speaker mentions that to run Kallisto for read alignment with paired end sequencing, one would enter the following:

where the sample files at the end refer to the two different end reads (forward and reverse read files). However, browsing the NCBI for paired end data, I noticed that only 1 Fastq file is provided per replicate for a paired end procedure (here's an example). So I'm a bit confused about this.

Building off of the answer by @An Ignorant Wanderer:
Essentially there are indeed the paired end data files available, but it takes a little bit of extra navigation to reach it.
Using the specific example you provided of BioProject PRJNA574273:
First select a specific sequencing "Run," aka the sequencing of a specific sample / experimental condition. For example, click on the link to run: SRR10189234
This brings you to the default tab: "Metadata." As @An Ignorant Wanderer pointed out, then you can navigate to the "Data Access" tab.
Then look under the "Original format" subheading. Here you'll find links to two files, corresponding to the two orientations of the paired end sequence files.
Edit: Also please see the answer by @burger, the website linked in that answer is much more user friendly in my opinion.

The problem is that the SRA website can be a little confusing. View the same project on ENA where the organization is more human-readable: https://www.ebi.ac.uk/ena/browser/view/PRJNA574273

Actually after clicking on one of the "Runs" and going to "data access" I noticed two files are provided per run

Paired ends is a configuration in sequencing platforms. Illumina is the most popular at the moment and most single-cell sequencing is done using paired-end Illumina sequencing. It varies between single-cell technologies but most 3’ scRNA-Seq protocols use paired ends.
Kallisto is a tool for handling raw data. Many databases supported processed data in BAM format or UMI/read counts. You should not need to run kallisto unless you are involved in handling raw data from sequencing performed by your laboratory. Some databases such as SRA and EGA support raw data in FASTQ format, others do not for ethical reasons or the data is too large.
I recommend downloading from SRA via their command line interface (this is documented on their website). It is possible to extract paired ends from an aligned BAM file using bedtools and samtools. I’ve done this with the Macosko et al (Cell 2015) DropSeq data for example.
Note that in some cases, the barcode is stored in databases as the headers of the FASTQ, BAM, BED, or CTSS file rather than a separate read. Single cell transcriptomics is an emerging field so there isn’t a standard way to do this, it varies between labs and data centres. For single-cell specifically, some databases are being set up, for example: https://data.humancellatlas.org/


text:
What unit (TPM, FPKM/RPKM, or other) to use when working across samples
I have raw gene read counts and would like to perform an analysis across multiple samples. I've found conflicting info online on how this should be done. One commonality however is that FPKM/RPKM should not be used. On the other hand, there seems (to me at least) to be conflict about whether TPM is suitable. For example, this StatQuest video seems to suggest that TPM is suitable for analysis across samples. One person on this Biostar post (among other resources I've looked at) suggested that none of the units (FPKM, RPKM, TPM) are suitable for cross sample comparisons.
I am quite confused. Can someone clear up the confusion?

Choices in bioinformatics tend to be quite situation dependent; could you please provide more information about the specific problem you're having? What tool are you using for mapping? Are the reads long reads, paired end, or single reads?


text:
Advantages of paired-end sequencing compared to single end
One of the advantages of paired end sequencing over single end is that it doubles the amount of data. Another supposed advantage is that it leads to more accurate reads because if say Read 1 (see picture below) maps to two different regions of the genome, Read 2 can be used to help determine which one of the two regions makes more sense. However, I don't understand how this second advantage pertains only to paired end sequencing. With single-end sequencing, couldn't you just use the read that comes right after read 1 to essentially get a longer net read?

An explanation for the downvote would be appreciated

You seem to think you can put the reads in order and they somehow join up, the reality is there is not order to the reads and it's only by assembling them like a jigsaw puzzle or aligning them to a reference genome that you can infer their location. I didn't downvote

@Chris_Rands thank you for your answer. Hmm, could you elaborate more as to why there is no order? I mean, when we're obtaining the reads, aren't we starting with a primer at one end of the fragment and working our way to the other end to generate the multiple reads that make up that fragment?

You have a picture of an Illumina read, which clearly does not work by 'working your way' from one end all the way through to the other.  Why technology are you talking about that works as you suggest, and why did you post a picture of an Illumina read when it's obviously not what you are talking about?  Why are you talking about multiple reads when you are asserting that there is an inherent ordering of single end reads?

There isn't any "order" that could be assumed from the FASTQ file.  That's why after you align reads, you often have to sort the generated SAM/BAM file.

@Ray so if there is no order to the FASTQ file, then how is it possible to infer which reads are paired together (in the case of paired end sequencing)?

Sorry for the confusion!  My "there is no order in the FASTQ file" comment refers to single-ended sequencing -- a single file.  For paired-end sequencing, the reads are usually in two separate files:  one for forward, the other for reverse.  Each read has a name, and the names are unique within a file.  Between files, they are matched up to form a pair by name.  So, yes, they are ordered in the sense that read called "ABC" is the 100th read in the forward file and it should be the 100th read in the reverse file.  A good program would also double-check the names match, too!

But within either file, there is no meaning between the 99th read and the 100th read (i.e., why it is the 100th read).  Thus, returning to single-ended sequencing (i.e., one file), there is no order that can be assumed.  I think some programs accept paired-end sequencing as a single file so that reads are interleaved.  Since they alternate, you can create two files from that easily, if you wanted to.

@Ray okay thanks this makes sense now

@AnIgnorantWanderer You're welcome!  Just to emphasize what someone else hinted, this refers to Illumina reads and how it (BaseSpace or equivalent programs) changes the images to FASTQ files.  Another sequencing technology might choose another format that is still FASTQ, but maybe order of reads exists.  (I don't know of one off hand, but I won't rule that out that it exists now or in the future.)  Good luck!

Like you mentioned, having read pairs can help with alignment.
Your FASTQ for a single-end sequencing run will look something like this:
@SEQ_ID_0001
AGCTAGCGCGGTTGGCTTAGCGACT
+
!''*((((***+))%%%++)(%%%%
@SEQ_ID_0002
AGGTGGTTGTAGGGAAAAAAGTCTC
+
!''*((((***+))%%%++)(%%%%
...

There is not necessarily a relationship between @SEQ_ID_0001 and @SEQ_ID_0002.
They could come from neigbouring regions in the genome or they could come from different species if that's what your experiment involves.
So aligning @SEQ_ID_0001 doesn't exactly give you information about where SEQ_ID_0002 goes.
@SEQ_ID_0001 could equally map to multiple locations in the genome with the same probability.
How do we decide where it "should" go?
We can randomly distribute multi-mapped reads or do some other fancy stuff, but that usually relies on facts about the data as a whole, and not something specific to the original DNA molecule that sequencing read came from.
If we have paired-end sequencing, then like your diagram shows, you know that read 1 (R1) and read 2 (R2; these have the same @SEQ_ID) will relate to each other in some way.
@SEQ_ID_0001 R1 may be multi-mapped, but if @SEQ_ID_0001 R2 isn't, then we know that R1 should probably be relatively close to R2.
Because of how you prepare your DNA samples before sequencing, you'll have some information about the insert length distribution.
You can use this information with sequence aligners to better identify where R1 should be placed.
You can also do this if both R1 and R2 are multi-mapped; certain pairwise alignments will be more likely than others, even if the marginal alignment probabilities for each read are the same.
Finally, there are other unexpected benefits that you can get from paired-end sequencing.
All of the above applies to just the mapping step, but applications like structural variant detection and chromosome conformation capture really benefit from paired-end sequencing.
If you have single-end reads for structural variant detection, then you will only detect breakpoints if the read overlaps a breakpoint.
But if you have paired-end, you not only have twice the amount of DNA sequenced to look for breakpoints, you can also estimate if there is a breakpoint between the two reads (in the inner distance from your diagram) if the two reads align really far apart from each other, or on separate chromosomes.
The idea is similar for chromosome conformation capture and the ligation junction sites.

In paired-end sequencing, you have an expected fragment length. So if Read 1 maps to two different places, whereas Read 2 only maps to one place. Then you can determine which one is more likely for Read 1 to map to. It cannot be too far away from Read 2 given the expected fragment length.
It is not necessarily true that paired-end sequencing gets twice the amount of data. The total amount of data can be determined when running sequencing. If single-end has 20M reads, I can just stop sequencing at 20M as well for paired-end.


text:
Difference between PPI networks and Gene Co-expression networks
What exactly is the difference between protein-protein interaction (PPI) networks and gene co-expression networks? Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks, and so represent processed data. But could someone elaborate more on the differences?

Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks

Not really, no. PPI stands for protein-protein interaction. In PPI networks, each node is a protein and each edge (line) connecting it to another node represents an interaction. A very common example is that an edge represents a physical interaction between these two proteins. This means that there have been experiments to indicate that these two proteins come into contact in the cell and interact.
For example, they might form part of a protein complex, so their interaction was seen using chromatin immunoprecipitation, where the two proteins were captured together. Or, they may have been seen together using FRET, a technique that lets you detect when two molecules are in very close proximity.
Alternatively, some interactions can be theoretical, where in silico modelling has suggested that these two proteins can interact, even though we don't have any direct evidence that they do so.
Finally, yes, it is also possible to make a PPI network based on gene co-expression data, under the assumption that if two genes are co-expressed, then their protein products will be present in the same tissue. However, that doesn't mean they interact in any way, so it isn't really a very common way of building PPI networks, as far as I know.
In co-expression networks, on the other hand, each node is a gene and not its protein product. And here, the edges indicate the coexpression. Often the level of correlation between the two genes' expression data.
The most important difference between these two kinds of network, apart from the obvious mentioned above, is that they are used in different ways. A gene coexpression network lets you identify genes that are regulated in similar ways. It is a tool you can use to investigate the regulatory network of the cell, to see what genes are "turned on" or "off" together.
The PPI network, on the other hand, is a way of understanding the result of this gene co-regulation. You are not looking at how the gene network is regulated any more, but how the products of the genes interact with one another. So it is a very different focus.
You can maybe think of this in terms of school children and their parents. You could make a network where each node represents a pair of parents, and then draw lines between them if their children were born in the same year. That would give you a general overview of when this group of people procreated. However, that two children were born on the same day doesn't necessarily imply that the children themselves will ever come into contact. This could be the equivalent of the gene co-expression network.
Then, make another network where each node is a child and draw a line between any children who know each other. This would be the equivalent of the PPI network. You would likely find some correlations between the two, with children who know each other being more likely to have been born at around the same time, but the two networks give you very different information.

Co-expression networks are based on gene expression data. Genes are compared using their expression profiles over multiple samples. A common distance measurement is the Pearson Correlation Coefficient, assuming that co-expressed genes (genes expressed simultaneously in similar conditions) will have a high correlation coefficient.
In its simple form, a co-expression network is done like this. An expression matrix containing n genes x p samples is used in input (from microarrays or RNA-seq) to calculate a distance matrix of size n x n. Then the distance matrix is thresholded to keep only the most significant relationship (eg correlation >0.8). Gene pairs are connected if their correlation coefficient is above this threshold. Then all connections are visualized in a network, which can then be subjected to community detection algorithms to find more densely connected set of genes.
PPI are usually built on in vitro assays such as Yeast-2-Hybrid. But connections can also be inferred from literature data and expression data. You may have a look at the StringDB for example to see which kind of data can be used to connect proteins.

It is also possible to infer regulatory networks from gene expression data using specific data (eg mutants in transcription factors and controls) and specific metrics (eg bayesian inference). Infering a regulatory network is generally more challenging.