text:
Identifying mutually **exclusive** gene sets
I am interested in identifying gene pairs (or better: sets of genes) whose expression is mutually exclusive.
Ideally, both genes (or gene sets) would be widely expressed but I am also interested in cases that are very cell type specific.
Familiar examples of the former would include some of the cyclin genes, in humans e.g. cyclin E and cyfclin B, which are expressed during S-phase and M-phase, respectively. Another example would be some of the clock-genes, where Per1 and Per2 are often highly expressed in the suprachiasmatic nucleus (SCN) before midday, and whereas Bmal1 is expressed after midday.
However, such genes need not be cyclically expressed, they can simply be very specific markers for certain cell types. For example, GAD65 and GAD67 are very specific for interneurons, whereas CamKIIa is very specific for excitatory glutamatergic cells.
I have been able to find a single paper so far that systematically tries to find such genes from gene expression data, namely: Tzanis & Vlahavas (2010) Mining for Mutually Exclusive Gene Expressions. According to Google Scholar, this paper is cited by a single other paper that does not seem to be relevant. In Tzanis & Vlahavas (2010), they use a SAGE library (from 2002) of short expression sequence tags to find these strong negative associations. By today's standards, the data set is pretty small (basically, it contains 90 samples). The paper itself is hence anything but a canonical source for such negative interactions.
Questions:

Is there a curated database of mutually exclusive gene pairs or gene sets?
If not, is there a curated database that would allow me to find such genes?

Otherwise, my current plan would be to use the human or mouse cell atlas single cell data to define a signed, weighted gene expression network.
I would then search for pairs of clusters/modules, where genes have predominantly positive correlations within a cluster but predominantly negative correlations between clusters.
So essentially, signed, weighted gene correlation network analysis (WGCNA) with extra steps.

Has somebody published such an analysis already (the WGCNA literature is pretty substantial)?
Is there other (recent) related work that I am unaware off?

I am fairly new to bioinformatics, so the lack of papers and resources that I am able to find is likely mostly a reflection of my lack of ability to find the correct search term.
Addendum
Similar biostars thread from 5 months ago, also without an answer.

"I am interested in identifying gene pairs (or better: sets of genes) whose expression is mutually exclusive." To what biological end? Are you looking for marker genes/sets? There's extensive literature on that.

@DevonRyan The most basic, direct application that comes to mind would be to identify doublets in single cell RNA sequencing data. Most existing approaches are data set imminent and hence unreliable when the number of samples is small. More broadly speaking, I am thinking about ways to bake in biological priors into RNASeq analysis pipelines. From my own literature review, mutually exclusive interactions seemed to me an under-explored direction of research in this regard, and hence a good starting point.

The only circumstance where doublets are an issue is when the technology is using microfluidics (or similar) and then you don't have issues with small sample sizes. It's unclear how mutual exclusivity (as much as that ever actually happens) would be relevant to generic RNA-seq analysis. As is this post needs much more focus on a single specific application so the solution can be appropriate, as there will exist no single generic solution.

@DevonRyan "The only circumstance where doublets are an issue is when the technology is using microfluidics (or similar) and then you don't have issues with small sample sizes. " -- Having access to (unpublished) data sets that are a) small, and b) may contain an appreciable number of doublets, I respectfully disagree. Irrespective of the sequencing method, not all tissues are easily dissociated such that you can get a single cell suspension that you could process further using conventional techniques.

@DevonRyan "As is this post needs much more focus on a single specific application so the solution can be appropriate, as there will exist no single generic solution." -- Why in your opinion would my questions, in particular question 1 & 2, require more focus? If I were to ask a different but analogous question, e.g. if there exists a database that lists relationships between signalling molecules or metabolites, everybody would direct me towards KEGG. I am not trying to be combative -- I am new on this forum and trying to understand what extra-information you would find critical.

The exact biological end goal is critical information. Just asking for databases in and of themselves is fine (FYI they won't exist), but once you start asking questions 2-4 you need information on exact context.

BTW, if dissociation is an issue and low-throughput methods are being used anyway then I'd suggest just using a microscope for dublet detection.

@DevonRyan We (have to) harvest the RNA in situ. We image the cells while we extract the RNA but that is not fool-proof.

I foresee a tie in with spatial single cell sequencing then. This is all relevant information that can be added to the post.

@DevonRyan The particular data set I am working on has been obtained using Patch-seq (basically). However, I foresee other applications, where you might have low throughput and high chance of contamination. For example, you might want isolate RNA from a synaptic boutton, and then you might worry about contamination of the synaptic (and hence neuronal) RNA from neighbouring astrocytes, which often engulf synapses with their processes. Since I think there will actually be quite a number of such applications, I tried to keep the question broad.

I do see your point that in particular Q4 does not have a canonical answer. I don't have time to rework the post now, but will do so soon.

"Mutually exclusive" is not a precise statistical term because it could be seen as independence.
Anyway if truly mutually exclusive AND their behaviour follows periodicity, i.e. sine waves, which it might do in a cell cycle, i.e. the cells divide at periodic intervals, you can do this via asynchronicity.
If the gene expression can be described under periodicity analysis you simply sum the amplitude between two genes. When the amplitude is zero ... they are expression is antogonistic - when one is expressed the other is silent.
Amplitude simply means the height of the period.
Periodicity analysis is quite simple if you have uniform time measurements.
If you have 'missing data', i.e. you time intervals are not uniform, for example you skipped a reading then this is called missing data in periodicity analysis. Missing data periodicity analysis is complicated, but its doable.
So in this context its quite simple, you would identify the periodicity of all genes and construct a 2x2 matrix subtracting their respective amplitudes. It would be kinda juicy to do.

The above is cool if the genes are periodic, however covariance is a better solution in general terms.
What you want is negative co-variance. If two genes had a covariance of -1 ... that will mutually exclusive. As above 2x2 matrix and the co-variance across the grid. You do need to be careful that the program you are using doesn't remove negative values of the covariance.
Covariance = 1, meaning two genes are co-expressed
Covariance = 0, random distribution
Covariance = -1 mutually exclusive expression


text:
When running ChIPQC bioconductor package after peak calling, should I use sorted BAM files or unsorted BAM files?
I have called peaks using MACS2 and now I'm trying to run ChIPQC bioconductor package. Can someone please tell me should I use sorted BAM files or unsorted BAM files? I have created sorted BAM files and from that indexes for BAM files. When I was trying run ChIPQC earlier, I got an error. Then I figured out that is because I don't have index files. Then I created indexes after sorting BAM files.

Please add relevant code and error messages. Anecdotal descriptions are never sufficient.

You need to provide sorted bam files, ideally with an index. If you run ChIPQCsample, it calls ChIPQC:::sampleQC which iterates through the chromosomes or contigs, using the index, and the absence of the index will crash it.
We can use an example, first in linux:
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibTfbs/wgEncodeHaibTfbsA549RxlchV0416101Etoh02AlnRep1.bam -O Etoh02AlnRep1.bam

samtools index ./Etoh02AlnRep1.bam
# s
samtools view -uh Etoh02AlnRep1.bam chr1 > sorted.bam
samtools index sorted.bam
# quick way to unsort
samtools sort -n -o unsorted.bam sorted.bam
# it fails to index
samtools index unsorted.bam
[E::hts_idx_push] Unsorted positions on sequence #1: 246814554 followed by 42922635
samtools index: failed to create index for "unsorted.bam"

Then in R:
library(ChIPQC)
ChIPQCsample("unsorted.bam")
Bam file has 25 contigs
Error in value[[3L]](cond) :
  failed to open BamFile: failed to load BAM index
  file: unsorted.bam

ChIPQCsample("sorted.bam",annotation=NULL)
Bam file has 25 contigs
Calculating coverage histogram for chr1

Calculating SSD for chr1

Calculating unique positions per strand for chr1

Calculating shift for chr1

 300 / 300
Counting reads in features for chr1

....

source: Stack Exchange

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

sort: votes


text:
Distinguishing Splicings that Produce Functional Transcripts Amongst Alternates Shown
This is a link to the UCSC Genome Browser showing a portion of Ch 17 containing TP53 here.
There are several representations of the gene of varying length and number of exons. I presume these are equivalent to alternate splicings.
Assuming this presumption is correct(?), how can I determine which result in RNAs that are experimentally verified?

@Michael G. Thanks, nice touch.

There are two resources worth checking: GTEx portal and Ensembl.
Uniprot and NCBI will call "canonical" the longest sequence comprising valid exons. These are not necessarily the biologically valid ones.
According to Ensembl (TP53, there are a bunch that are experimentally validated and are all class 1. Note the Uniprot column being nearly homologeous: to match the transcripts one has to look at the length, in this case TP53-201 is 393aa long, which is class 1. As mentioned this is not always the case.
So which are more expressed in tissue X?
In GTEx, under the collapsed div labelled "Exon Expression" press the "isoform Expression" tab. The most expressed is ENST00000620739.4 (354aa long) join with ENST00000269305.8 (393aa).
In the Exon Expression you will find out what exons are the issue.
Now. to get the sequence of a transcript or covert a variant Ensembl drops the ball. There is no clear web tab that I know of. Nor an API route.
ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/pep/Homo_sapiens.GRCh38.pep.all.fa.gz contains all the sequences however. For a project of mine (ENSTMapper in this link) to map an amino acid mutation on one transcript to another I ended up write my own short converter, although there are a few similar tools out there.

(Were this answer to be satisfactory, no need to mark it as such: it is not really worth 50 points —keep them as they unlock useful privileges. An upvote is fine —and is way more than most users)

Many thanks for your generosity. But you put a substantial effort, and I like to stick to my word. Best regards,

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


text:
Reading multiple raw files in Seurat
I have multiple single cell samples to analyze and I'm following the instructions in Satija Lab's website. I want to merge all the count files from all the samples at once, and associate the metadata to each sample.
For the Read10x command, is there a way to read multiple files at once, such as merging all the filter_featire_bc_matrix files?
I have all the samples saved in a directory on my desktop, is it possible to call the directory within the script for it to read all the raw counts?
This is what I have thus far:
# load the required libraries
library(dplyr)
library(Seurat)
library(ggplot2)

-----------------------------------------------

# set up the working directory

wd = ("~/Seurat") setwd(wd)

# read all metadata

metadata = read.table("metadata_Jan2020.csv", sep=';', header=T, stringsAsFactors=F)
rownames(metadata) = metadata$LibraryID

#metadata
#mf = metadata

mf = metadata %>% filter(Type %in% c("PDAC", "NORMAL"))
mf = mf %>% filter(!(ID %in% c("PDAC-HTB2854")))`

Seurat has a vignette specifically for combining multiple 10x libraries.

Yes, I know. But I'm wondering if there is a way to merge all counts before creating a Seurat object. I'd like to just create one Seurat object with all the counts.

You don't have to keep the intermediate Seurat objects.

To merge all counts before creating individual Seurat objects, you will need to give a prefix or a suffix to cell names. I would however advise to create individual Seurat objects with apply() or mclapply() and then reduce() these with Seurat's merge(), this will give you a single Seurat object with all your samples.

@haci How do I give a prefix or suffix?

Note that a month later, they put out a new vignette about merging data sets which is more complex:  https://satijalab.org/seurat/pancreas_integration_label_transfer.html

That vignette is over a year old. Also, it is for integration (batch correction), not simple merging of data sets.

That's not what I'm looking for. Is there not a way to merge the counts from a directory and then create just one Seurat object?

See previous comment by @haci. However, there are many benefits of creating separate Seurat objects for each library, such as performing per-library QC.


text:
What are the different kinds of bioluminescent genes?
I know of the common green glow gene but I forgot the name and I also know that some algae glow blue. There are so many types of bioluminesent organisms, so I am wondering what species have which genes and what other genes are associated with their metabolism and the systems that regulate these genes.

Hi, I think this is better moved to se.biology

I was hoping to compile a list of answers with genes and related species, or I could comment about the species after searching the BLAST database for each gene. Maybe someone would answer with a sequence that I could look up on the BLAST database. I need a review of discussion about genes and bioinformatics, and thought I would start with bioluminesence. Since some time has elapsed since I studied this, I was thinking more information has been developed about this topic, particularly Luciferin and related species. Answers shedding light on metabolism-related studies of genes are desired.

BLAST is a search algorithm, do you mean NCBI database? Would it make sense to reformulate the question as "how to find species with bioluminiscence and preferantially also parthways involved in the process?" - because that would have a good potential to have a decent bioinformatics answer. Note you can [edit] your question (and it's better than writing long comments)

There are three processes wherein light is emitted.

Bioluminescence: a chemical reaction releases light. The enzyme that does this is luciferase, while its substrate is luciferin, a small molecule. There are half a dozen organisms that have analogous (evolved independently) enzymes and totally different substrates, such as fireflies, Alivibrio bacteria, diatomes, deep sea fish, some mushrooms etc. The wiki article has a list of these and their small molecules.
Fluorescence: light is absorbed by the molecule and is emitted at a lower energy (wider wavelength, redshifted). Green fluorescent protein of Aequorea victoria jellyfish is the classic example of this. It is excited with UV and emits green light. Here is an interactive model, I threw together. Variants exists for biotech reasons that span the whole of the visible spectrum. Other species have a homologous beta barrel, such as corals (red).
Phosphorescence: light is absorbed by the molecule and emitted later.  Ignoring the picosecond delay in fluorescent small molecules (cf. fluorescent lifetime microscopy), this is exclusively a properties of lanthanide elements (e.g. Europium). Lathanide binding domains have been designed, but do no exist in nature.

Thank you for classifying types of bioluminesence and connecting some species. I am wondering what the genes are so that I can search for the species. For example, many jellyfish, deep sea species like the lanternfish, and algae glow, and the cuddlefish utilizes chromatophores to change color. I dont think you mentioned fungus, any idea what kind of glow gene they have?

Cuddlefish chromatophores are simply pigments, i.e. small molecules with a resonant structure but do not release the excited energy as a photon.

Fungal bioluminescence was described fairly recently (https://advances.sciencemag.org/content/3/4/e1602847.full). There are two parts the enzyme and the pathway to synthesise the substrate (3-hydroxy hispidin). The pathway biosynthesis is quite clear, but https://pubs.rsc.org/en/content/articlelanding/2010/np/c004593d/unauth#!divAbstract discusses some enzymes. N. gardneri is not sequenced, but other funghi are.

In terms of gene sequences. You'll just have to do a search for luciferin, which is a generic name. Although most may not be sequenced, eg. glow-worm genome.

Right, I am figuring that out. How many species have been genomically evaluated, and just specifcially because they glow? rhetorical question) Especially because derivative molecules have already been innovated in labs, there isn't really a point. I think I was hoping that since it has been a few years and gene sequencing machines are easier to access, that more species would be on file with more data about their molecules and metabolism of the genes that regulate them.

Well, it is not trivial finding the enzymes that encode a secondary metabolism pathway, so sequencing is not the discouraging step. In thermophilic bacteria they will be in the same operon and in fungi there will be some clustering, but elsewhere its horrid. Scientific effort goes into seeing what the small molecule substrate is and if it could be technologically useful. If it is, such as firefly (very bright) or bacterial (easily usable pathway in SynBio), it will get more attention.

The green protein is called GFP of 238 amino acid residues (26.9 kDa) and produced by:

the jellyfish Aequorea victoria, It has a major excitation peak at a wavelength of 395 nm.
the sea pansy, Renilla reniformis of 498 nm.

However, the protein is widespread and not restricted to these species.

Green Florescent Protein, thanks. Do you know the 'Ascension number' of the gene or where I can find the sequence?

Looking up the sequence, I found this site: https://www.uniprot.org/uniprot/P42212 that had  one listing for this gene here: https://www.ebi.ac.uk/ena/data/view/M62654 so then I looked up "M62654" in the NIH's BLAST database to pull up a list of 100 results:   https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome    There were 5 from the jellyfish, one from a species Azotobacter vinelandii, and a listing for a synthetically derived version from E. coli. There were several related genes that were vectors, like a promotor vector at the end of the list

There are also other versions of GFP (such as "enhanced GFP") that are more commonly used in the lab.

There appears to be a dinoflagellate with Luciferase and the diversity of their genes was studied in the Baltic, found described in an article here: https://www.tandfonline.com/doi/full/10.1080/09670262.2016.1160441
This article by Smithsonian: https://ocean.si.edu/ocean-life/fish/bioluminescence describes a lot of Luciferin metabolizing species in the oceans, saying that different species independently evolved the forms they use. It describes some of the systems like endosymbiotic bacteria, special tissue, and photo-proteins.


text:
What does this statement mean?
I have a copy number segment file with this by SCAT R package
Chromosome  Start   End
1           10583   120983
1           121009  245455
1           404187  5719467
1           5719471 5735807

I want to count the number of markers in these locations
I asked somebody how to do this, he replied:
We created marker file by dividing the genome into 10 kb bins, so the markers file contains the start and stop positions for each of these 10kb segments for each chromosome along with ASCAT segments start and end positions. Then, for each ASCAT copy number segment, we count how many markers fall within that region for total probes.
Can somebody helps to understand what he mean by dividing the genome

Why don't you ask the individual who gave that advice instead of asking others to interpret it?

Thank you, I emailed him but he did not reply me back

They created 10kb bins, which could be called "dividing the genome". These are just 10kb long contiguous stretches (so position 1-10000,then 10001-20000, etc.). This is, of course, a nonsensical way of going about things, because "counting those up" in the regions you posted is the same as taking the width of the regions, dividing by 10000, and rounding a bit.
Presumably you have your markers in something like BED format, so just use bedtools to intersect them and do the counting for you.


text:
Looking for a tool to find 16S rRNA in hundreds of genomes
I have hundreds of bacterial genomes and I would like to know which is the best informatic tool in which I can run all genomes on one go and get a 16S rRNA gene sequence of each genome as output.
Could anyone give me some advice?

Have you tried anything any particular? As a general rule, there is no single "best" tool for much of anything, you just need something to give you a reliable result.

You can use Barrnap
barrnap -k bac --threads "$(nproc)" \
        -o output_rrna.fna < input_file.fna \
        > output_rrna.gff 2> /dev/null

Put the file in one line of sequence per fasta id and get sequence with 16S annotation.
perl -pe 'if(/\>/){s/\n/\t/}; s/\n//; s/\>/\n\>/' output_rrna.fa | \
    perl -pe 's/\t/\n/' | tail -n+2 | cut -d\: -f1 | \
    grep -A1 16S > output_16S.fasta

Check the output_rrna.gff and make a BLAST of the sequence to be sure about the predicted 16S.

Thanks zorbax, I am going to have a look and I will give you a comment about what I get

You can perform genome annotation using tools like prokka or PGAP. Using some python scripts  you can extract the predicted 16S and blast to confirm.


text:
Can someone recommend software for reference-based genotype-estimation?
Problem: I have low-coverage sequencing data from humans in vcf format. For each sample, I have PL scores / preliminary genotype calls from the software ATLAS (it's a bit like GATK for those who don't know). Because the samples are low coverage (mean 1x, with lots of missing data), I am looking to enhance the genotype calls by using a reference panel and genetic map, before go on to phase the samples. Briefly, the reference panel and genetic map can be used as a kind of prior to inform one about genotypes when there is little data present.
What I've done so far: So far I have been using Beagle4 with a command:
beagle4 gl=lowCoverageSamples.vcf ref=highCoverageReference.vcf map=geneticMap out=lowCoverageSamples.newCalls.vcf

However, Beagle4 is extremely slow, very memory inefficient and quite old. There have been many advances in phasing methods, which are now much more accurate, faster and memory efficient than the beagle4 model. However, I haven't come across any similar improved methods for the genotype re-estimation step.
Does anyone knew of any newer pieces of software which do such a thing, which may more faster and more accurate than Beagle4. It needs to have 2 things in particular:
1) Accept data input in PL score format, as to allow the appropriate representation of uncertainty in low coverage samples.
2) Produce a genotype call / genotype probability for each SNP in each individual as output.

In case anyone was interested, a new piece of software for rapidly imputing and phasing low coverage sequence data from genotype likelihoods has just been released. It is called GLIMSPE and is exactly what I was after.
It looks to be much faster (1200x speedup) compared to Beagle 4.0, and is more accurate at imputing low frequency variants. It is also more amenable to use with larger reference panels (such as the Haplotype Reference Consortium panel).


text:
Database with antigens and current vaccines
I am seeking virological databases focussing on:

vaccines
antibody interactions
host receptors.

Preferably the data base will host multiple strains/samples of a given virus samples.
Any assistance would be welcome.

Hi MniNi, could you tell us what viruses you are interested in?

Anything that is parse-able data, has documented binding site or vaccine, best both. I need real data in some data analysis I want to try on, with sequential substituted correlation. (Having substituents on certain shapes of chain, and correlating differences to differences in those substituents)

Hey ya, note you can [edit] the question and add more details (as those you wrote in the comment).

Given "any virus will do", the viruses that are "head and shoulders" above the rest in terms of volumes of data are HIV and influenza.
The influenza database is here
The influenza viral spike that attaches to the cell receptor is the HA protein – hemagglutinin. The cell receptor is sialic acid – a small sugar that is attached to many different proteins on the cell surface. You can of course simply search PDB for HA or else Genbank.

Ouch, I was abuncted by the goverment again, so I think somebody took the best from influenza vaccine technique and applied it to corona like I wanted.


text:
vcf file to gene annotations and or genes
Is there some tool that will convert a vcf/bcf file to some gene annotation, based on the fact that I have the gene annotation of the already exsisting reference genome. Or some tool that will give me the genes based on the most probable chance.

VCF files contain information on variants. Each position in a VCF file can be attributed to a gene (if it falls within a gene), is that what you're looking for? By "gene annotation of the already existing reference genome", are you referring to a GTF/GFF file?

I believe the question is about using VEP. https://www.ensembl.org/info/docs/tools/vep/index.html maybe?

People should really stop abusing the word "convert". Does OP wish to annotate their VCF or query a larger dataset based on some data in the VCF or do they think they can actually convert information about variants to information about genes?

I already have annotate my VCF(actually bcf file) file with SnpEff. So I sort of doubt whether it was possible to convert it to information about genes. On one hand it would be logically if it could on a biological perspective, on a technical perspective however, it wouldn't be that strange it that isn't that easy.

Yes I think I need VEP

VEP and snpEff give you similar information, as does ANNOVAR. What is the "information about genes" that you are looking for?

Like what the actual sequence is mutated and by reference

Creating a FASTA file by applying variants found in a VCF file is called "creating a consensus FASTA file". You can use bcftools consensus. Unfortunately, the tools only works genome-wide (at the FASTA level), so you might need to use a different tool to extract a gene sequence from the FASTA file using the gene's co-ordinates.
If the VCF file is annotated with HGVS notations for the variants, there might be a way to use those notations to apply the edit to the reference gene/transcript sequence directly. I know the hgvs python module works well with the HGVS notation, but I don't see a way to "apply" the variant to the reference sequence.


text:
How to download multiple proteomes at once?
I'm looking for a way to download multiple proteomes at once from one clade, as fasta. Probably from the NCBI, because it looks more user-friendly than Ensembl, but Ensembl is okay too. In the best case scenario, I would like to constrain the download to a specific clade (eg. Actinobacteria, Chloroflexi etc.), maybe with a specific RefSeq Category (eg. reference/representative genome) or even depending on the genome assembly level. Command line would probably be the best option, but if there's a direct download option it would be nice to know it as well.
Does anything like this exist? Help would be much appreciated!

The FTP has it: ftp://ftp.ncbi.nlm.nih.gov/genomes/ therein look for the fa.gz file

You could try download-refseq-genomes. It will fetch all genomes from the NCBI FTP server that are in a specific subtree of the phylogeny.
For example, downloading the amino acid sequences from all proteins in each genome assembly for all species belonging to Chloroflexi (NCBI taxon ID = 200795):
download-refseq-genomes.pl -t faa 200795

This will download all assemblies with status=Complete. Use option -a to download all assemblies with any status.

You can use NCBI resources to download these data directly. For the entire proteome, you can download the sequences of all proteins annotated on the genome assembly in FASTA format. In the NCBI FTP directories, these files have the names with the suffix _protein.faa.gz. More details about the files in the FTP directories are described in the README.
Command-line method using EDirect
This method uses Entrez Direct to first query the Assembly database for all available genomes, parses the output to extract the ftp path and downloads data using wget command.
esearch -db assembly -query 'txid200795[organism] AND lastest_refseq[filter]' \
  | esummary \
  | xtract -pattern DocumentSummary -element FtpPath_RefSeq \
  | while read -r url ; do
    path=$(echo $url | perl -pe 's/(GC[FA]_\d+.*)/\1\/\1_protein.faa.gz/g') ;
    wget -q --show-progress "$path" -P protein_data ;
  done

Here, I am restricting the results to latest RefSeq data by using the latest_refseq[filter] but you can use any other Entrez filters to further narrow down your results. Along the same lines, you can modify the code minimally as shown below to download multiple file types:
esearch -db assembly -query 'txid200795[organism] AND latest_refseq[filter]' \
  | esummary \
  | xtract -pattern DocumentSummary -element FtpPath_RefSeq \
  | while read -r url ; do
    base_path=$(echo $url | perl -pe 's/(GCF_\d+.*)/\1\/\1_/g') ;
    for suf in genomic.fna.gz protein.faa.gz ; do
      wget -q --show-progress "$base_path$suf" -P ncbi_data ;
    done ;
  done

GUI method
You can use the same query as shown above in the NCBI Assembly portal and then use the 'Download Assemblies' button to pick and choose the type of files you would like to download. By default, if you do not select any assemblies in the results, this interface downloads files for all results.


text:
Importance of Proper Pairs vs Aligned Reads for RNASeq data
I have stranded, paired end RNASeq reads that I have aligned using STAR. I plan do conduct a differential expression analysis with DESeq2. After running quality control checks, a good portion of my reads are aligned but I'm a little concerned about the proportion of properly paired reads

Here is the inner distance plot for the samples

For human RNA sequencing, does anyone have a rule of thumb that they'd consider is an appropriate proportion of properly paired reads? Are there any good papers that talk about this issue? As MultiQC marks the column as yellow, I'm concerned the proper pair percentage may not be satisfactory.
How important are properly paired reads for differential expression analyses?

Differential gene analysis sees only the read counts. So the proportion of not properly paired reads will have an effect on the counts, if the programme / software used, counts fragments / read pairs instead of single reads. For example HTSeq will count

From https://htseq.readthedocs.io/en/release_0.11.1/count.html
For paired-end data, does htseq-count count reads or read pairs? Read
  pairs. The script is designed to count “units of evidence” for gene
  expression. If both mates map to the same gene, this still only shows
  that one cDNA fragment originated from that gene. Hence, it should be
  counted only once.

In Rsubread, featureCounts you set requireBothEndsMapped=TRUE to count fragments.
Short anaswer is if the sequencing is deep enough, you should still have enough reads. In your situation, I estimated that you lose about 40% of your reads if you use only properly paired reads, which means you lose quite a bit of counts.
The larger question is what reads are you missing?
I see that you have about 20-30% of reads that are non-uniquely mapped, I suppose part of this will contribute towards the not properly paired reads. I would check first whether most of the not properly paired comes from multiple aligned.
Then if not, you can check where these not properly paired reads are falling, if they are some repetitive regions or say rRNA etc, then most likely thats a problem with the library prep..

If you check the RSeQC docs you can get an explanation of how the inner distance stats work. By the look of your inner-distance MultiQC plot I guess that your sequencing is 2x100bp reads. This plot shows that your insert sizes on the pairs is pretty short, meaning that a lot of the paired reads are overlapping significantly. This will be affecting the alignments a little, resulting in the low proportion of properly paired reads.
The good news is that for differential expression, this probably doesn't really matter too much. What is important is getting fragment counts, so as long as at least one of the pair mates is mapping properly, then you will will be getting sensible counts - even if read 2 is not mapping properly. You should be able to have confidence in your results from multiple replicates if in doubt - you have more than just -rep1 right? ;)
Note that the yellow background in the MultiQC column is not passing judgement with categories like pass/warn/fail in FastQC - it just happens that the colour scheme ended up there for your samples. It's a continuous gradient from red to yellow to green for those values (0-100%).


text:
Trying to divide number in fourth TSV column by a certain number according to its corresponding identifier in the first column
I was wondering if anyone could help with my little issue in manipulating tsv files in order to calculate estimated coverage across many tsv files.
I have tsv files that look something like this:
ID        Info    Info2  x
M80428    info    info   1000
T83241    info    info   8400
E81292    info    info   10200

I also have another tsv file with all possible IDs and a corresponding number like so:
ID        y
M80428    27311
E11123    46531
P38163    14213
A34197    12892
T83241    37416
Q21345    27321
E81292    13101

x in the first tsv are estimated read counts and the y in the second is the genome length of each ID.
I would like to divide the "x" by the "y" in the second tsv that has the corresponding ID in order to get the estimated coverage of each hit. e.g. the first row would be 1000/27311, second row 8400/37416. I would then want this to be written to the fifth column of the tsv.
I am working in bash shell mostly here, and think that the way I would probably do this is with awk and by making the second tsv an associative array but I can't work out how exactly to do so. The python package dict seems like another option but my Python skills are practically zero.
Thanks for the help in advance!

I'm voting to close this question as off-topic because there is nothing related to bioinformatics in this post

Hi @stanleyHo, could you provide more detail, what have you tried and what biological context is your problem?

x in the first tsv are estimated read counts and the y in the second is the genome length of each ID. Trying to divide them to get estimated coverage.
I have tried to set the second tsv as variables and then tried to awk out the columns with divisions but i'm stuck at how to link the ID and y together.
Something like this:
cat genome-size.txt | while read eachline;
do echo $eachline;
vname=echo $eachline| awk '{print $1}';
div=echo $eachline| awk '{print $2}'
cat .tsv | awk -v div="$div" '{FS ="\t"} {OFS ="\t"} {print $1,$2,$3,$4,$4/div}' > new.tsv;
done

Hi @StanleyHo, note you can [edit] your questions to make them clearer.

Use Python's pandas ...

Python dictionaries are indeed a good tool for that.
Here is a script that takes 2 arguments (the file containing the counts and the file containing the genome length info), and prints the input data, with an added column containing coverage:
#!/usr/bin/env python3
import sys

counts_fname = sys.argv[1]
genome_lengths_fname = sys.argv[2]

genome_length_dict = {}
with open(genome_lengths_fname) as genome_lengths_f:
    # Checking the header line, for extra safety
    header = genome_lengths_f.readline().strip().split("\t")
    assert header[0] == "ID"
    assert header[1] == "y"
    # Reading the data lines
    for line in genome_lengths_f:
        [genome_id, genome_length] = line.strip().split("\t")
        genome_length_dict[genome_id] = int(genome_length)

with open(counts_fname) as counts_f:
    print("ID\tInfo\tInfo2\tx\tx/y")
    # Checking the header line, for extra safety
    header = counts_f.readline().strip().split("\t")
    assert header[0] == "ID"
    assert header[1] == "Info"
    assert header[2] == "Info2"
    assert header[3] == "x"
    # Reading the data lines
    for line in counts_f:
        [genome_id, info1, info2, counts] = line.strip().split("\t")
        # /!\ It would be different in python 2:
        # different print and division behaviours
        print(
            genome_id, info1, info2, int(counts),
            int(counts) / genome_length_dict[genome_id],
            sep="\t")

sys.exit(0)

Testing it:
$ ./coverage.py x.tsv y.tsv
ID  Info    Info2   x   x/y
M80428  info    info    1000    0.0366152832192157
T83241  info    info    8400    0.22450288646568314
E81292  info    info    10200   0.7785665216395695

This worked perfectly thanks so much!

I would do it in R:
tab1 <- read.table('file1.tsv', header = T, stringsAsFactors = F)
tab2 <- read.table('file2.tsv', header = T, stringsAsFactors = F)
merged_tav <- merge(reads_tab, genome_len_tab)
merged_tav$z <- merged_tav$x / merged_tav$y
merged_tav

returns
      ID Info Info2     x     y          z
1 E81292 info  info 10200 13101 0.77856652
2 M80428 info  info  1000 27311 0.03661528
3 T83241 info  info  8400 37416 0.22450289

note that only IDs that are in both tables will be preserved.


text:
Program Quality control for data sequencing (nanopore)
Hope all of you are fine. I am wondering if exists a program easy to use that can bring me Basic Statistics of Quality Control from raw data made by MinION device (nanopore).
I have read about FASTAQC and Porechop. My sequencing raw data was a result of a metagenomics analysis (16s and 18s). My main objective is only to filter the sequences and obtain a basic statistic about the lenght and the number of sequences ordered by the lenght. Thank you ladies and gentlemen
PD Sorry for my poor english. Day to day i'm trying to improve it.

A simple start would be to use MinIONQC, which is also available via bioconda under the name r-minionqc.

I assume raw data is in Fastq format, but ideally you should also use the sequencing_summary.txt file generated by the basecaller. Other alternatives for QC are PycoQC and NanoPlot. I developed and maintain NanoPlot, so let me know if there is something I can help you with.


text:
efetch nucleotide -> protein ids
I'm downloading translated viral genomes initially via Blast - which had shortcomings - and now by efetch.
What I want to do is obtain a complete list of all translated protein ids from a nucleotide Genbank entry.
efetch -db nucleotide -format gpc -id MT039887 | xtract -insd ....

The protein id and the translated gene sequences are within the efetch -db nucleotide -format gpc -id MT039887 output ... but how do I use xtract to parse them?

If you want the cds in fasta format, try:
esearch -db nucleotide -query MT039887 | efetch -format fasta_cds_na

And for the proteins sequence in fasta format:
elink -db nucleotide -target protein -id MT039887 | efetch -format fasta

Wow!! and thank you. Thats the answer ... days like this BioSE works. @zorbax I owe you. "Acknowledgements" thank you heading your way (thats not a metaphor).

In case you are still interested in just getting the protein ids and not sequence, you can use xtract like this:
efetch -db nucleotide -format gpc -id MT039887 | xtract  -insd CDS gene product protein_id
MT039887.1  orf1ab  orf1ab polyprotein           QHZ00388.1
MT039887.1  S       surface glycoprotein         QHZ00389.1
MT039887.1  orf3a   orf3a protein                QHZ00390.1
MT039887.1  E       envelope protein             QHZ00391.1
MT039887.1  M       membrane glycoprotein        QHZ00392.1
MT039887.1  orf6    orf6 protein                 QHZ00393.1
MT039887.1  orf7a   orf7a protein                QHZ00394.1
MT039887.1  orf8    orf8 protein                 QHZ00395.1
MT039887.1  N       nucleocapsid phosphoprotein  QHZ00396.1
MT039887.1  orf10   orf10 protein                QHZ00397.1

Thank you, this is really useful.


text:
Do additional peaks in percent GC of PacBio gDNA reads indicate contamination?
I have two sets of PacBio reads from genomic DNA of an Aspergillus species that were made from separate preps of  the culture. One of them has two additional peaks at 38% and 60% in the percent GC histogram produced by FastQC. Do these additional peaks indicate that there is contamination in that sample?

Yeah the prep looks odd, but you cannot really tell its contamination. You have to try and map the reads first, and see how many of them are unmapped, and see where the unmapped reads come from. You can just try with 1000 reads, should be obvious if it is contamination

Remake the graphs using only mapped reads. Then you can decide on your next steps.

The logic I would us is that if they come from the same genome, coverage peaks will go in multiplications of 2 - 1n, 2n, 4n, ..., while if it's a contamination, the coverage will go wild. So, I would recommend to look at the coverage vs GC to decide where it comes from.
You could either assemble the genome and use something like Blobtools to check coverage, GC and putative origin of sequences (and you should do that anyway).
Perhaps a bit simpler would be to look at it using k-mers, but that assumes you have reads with low error rate (In your case PacBio CCS or HiFi). If you do, you can use KAT that can do GC vs coverage plot - they intended it to look at GC bias in sequencing, but the vary same plot should show you if the three GC peaks have the same coverage or not.

I would suggest plugging your reads into a fast k-mer matching tool like mash screen to get a basic idea of where your reads come from, if you suspect contamination.
An issue with just looking at GC content is that some organelles/plasmids/etc. such as mitochondria have quite divergent GC content, and the different populations you observe might just indicate random fluctuations of copy number of these DNAs due to different sample preps. Not sure that would explain your plot but it does give me pause.
As suggested above, assembly is obviously a good idea if you are ok with the compute cost, and then you can use something like SPAdes or blobtools to pull out contigs from different origins.
KAT and similar tools as suggested are handy for hypothesis free evaluation of different populations in your reads, but they don't directly address the issue of contamination.
And obviously measure the quality of the reads to establish that they are ok and nothing weird happened there, they aren't all Gs or something.

Yes, it looks odd, but you might need to assemble the sequences to tell what is happening. The fact that the peaks are symmetrical around the presumably Aspergillus peak around 50% GC is strange. Usually contamination appears as one peak or shoulder on either side of the host DNA peak rather than symmetrical peaks.


text:
How to deal with zeros in the denominator when calculating percent enrichment?
I have a data.frame with results from a CRISPR screen.
I want to calculate the % enrichment of a variable of interest var from baseline to treatment group.
Here is the distribution of var in linear:

I can do scale var with log transformation

But many rows with var = 0 will be lost.
Moreover, I cannot compare % enrichment by dividing treatmentgroup::var to Control/baseline::var because I will run into the problem of dividing by zero.
I've seen people circumvent this problem by adding a small value ex)0.01 to the entire column of var.
But I'm wondering if there's a better way of circumventing this problem. For example, a simple transformation that can deal with zeros and negative values.
Appreciate any input!

The basic premise that I understand is this: There should not be a zero in your denominator. Zeros in the numerator can be circumvented by pseudocounts, but zero in the denominator indicates some sort of problem in the logic. Also, why are you working on CPMs directly? Why not use a differential expression tool?

@RamRS
Many thanks for the input - I'm a first year graduate student and have never worked with crispr screens. Right now I'm kind of brute forcing and see what I can get. Have you personally worked with any differential expression tool with one that you recommend trying out? And is MAGeCK one of those tools? Thanks so much again, appreciate any guidance!

I have worked with DESeq2 but I'm not sure if and how it would apply to CRISPR. Does current research in the field use count data for comparison? How does it account for library depth and other experimental factors?

@Nomad420 you may want to give a try to MAGeCK as you have mentioned. It accepts count data as input and will make your life easier.

Cube root transformation is another way to deal with zeros without dropping them.
var^(1/3)

EDIT as Michael pointed out, square root transformation is another--perhaps more common--transformation that leaves zeros as zeros.
sqrt(var)

Whilst I don't disagree, normally a square root transformation is used


text:
Joining columns to a sorted file
I have had a segmentation file (copy number) like below
  Sample      Chromosome  Start      End          Segment_Mean
LP6008334-DNA_C02   1     10583      62271      -0.992452416
LP6008334-DNA_C02   1   17201069    17230953    -0.992452416
LP6008334-DNA_C02   1   17235034    17235196     0.592510084
LP6008334-DNA_C02   1   142781738   142967188   -0.992452416

I was using bedtools to count the number of markers in above genomic range (Start, End) when I got as below
1   10583   62271   6
1   10583   62271   6
1   10583   750153  75
1   14907   1062025 106
1   121009  245455  13
1   404187  5719467 532
1   732809  5726842 500
1   821926  1120377 31
1   952428  1179983 23
1   1181751 1469134 29
1   3458681 143542402   14010
1   5726844 5736229 2
1   5735836 16832446    1111
1   5736682 8081290 236
1   8353020 8476455 13

Basically the same files but the second one has extra last column for the number of marker but sorted one I then I can not matched this with the above sample names and Segment_Mean
How I can add the corresponding sample names and Segment_Mean to my second file?
I have tried this where I am getting this mess
DF <- left_join(df_res1,df, by = "Start")

> head(DF)
           Sample.x Chromosome.x     Start     End.x Num_Probes          Sample.y Chromosome.y     End.y
1 LP6008334-DNA_C02            1     10583     62271          9 LP6008334-DNA_C02            1     62271
2 LP6008334-DNA_C02            1     10583     62271          9 LP6008337-DNA_A07            1     62271
3 LP6008334-DNA_C02            1     10583     62271          9 LP6008460-DNA_F02            1    750153
4 LP6008334-DNA_C02            1  17201069  17230953         20 LP6008334-DNA_C02            1  17230953
5 LP6008334-DNA_C02            1  17235034  17235196          1 LP6008334-DNA_C02            1  17235196
6 LP6008334-DNA_C02            1 142781738 142967188         40 LP6008334-DNA_C02            1 142967188
>

Also with this I have some overlaps like
300 segment overlaps detected in file '/temp/hgig/fi1d18/man/examplefiles/res1.txt'.

First overlap detected between segments at lines 1 and 9.

What have you tried? This is a problem you have faced and solved a few times before.

Thank you, I have edited main post with what I HAVE TRIED

Rearrange your text file into a sorted BED file:
$ tail -n+2 segmentation.txt | awk -vFS="\t" -vOFS="\t" '{ print $2, $3, $4, $1 }' | sort-bed - > segmentation.bed

Use BEDOPS bedmap with --echo-map-id-uniq and --count operators on merged regions:
$ bedmap --delim '\t' --echo --count --echo-map-id-uniq --fraction-map 1 <(bedops --merge segmentation.bed) segmentation.bed > answer.bed

The file answer.bed will contain merged intervals from segmentation.bed in the first three columns, the number of intervals from segmentation.bed which are contained entirely within the merged interval in the fourth column, and a unique listing of IDs of those contained intervals in the fifth column.

Sorry segmentation.txt is my output from DF <- left_join(df_res1,df, by = "Start")  ?

Sorry, you said you were using command-line tools, so I assumed that was the case here. Can you use write.table or similar to write data out of R to a text file?

I have already written DF dataframe drown by DF <- left_join(df_res1,df, by = "Start") as a txt file and uploaded that in where terminal can find that.

Sorry says [fi1d18@cyan01 examplefiles]$ tail -n+2 segmentation.txt | awk -vFS="\t" -vOFS="\t" '{ print $2, $3, $4, $1 }' | sort-bed - > segmentation.bed
Non-numeric start coordinate.  See line 312 in -.
(remember that chromosome names should not contain spaces.)
[fi1d18@cyan01 examplefiles]$    This is the link of my whole DF file  https://www.dropbox.com/s/54923xx5ds6p52d/DF.txt?dl=0

Maybe start with exporting data before the left join step, i.e., the data you posted at the top of your question? While I can make an informed guess, I'm not 100% sure what left_join does in this case, and I probably don't have time to dig into R and debug its syntactic idiosyncrasies to solve this there. I do think the command-line approach that I suggested will help you solve your problem, as described, however. Good luck and I hope this helps!


text:
Why does samtools mpileup sometimes include ref bases (other than ',' or '.')?
This is my first post here. I can think of no way of giving an easily reproducible example, as per the stack ethos. SO apologies in advance and any feedback on question format appreciated.
I have called samtools mpileup (Version: 1.9 (using htslib 1.9)) on some oxford nanopore data (mapped to GRCh38).
I am seeing some instances like below (grouped with pythons Counter), here reference is T at this position:
 Counter({',': 34, '.': 33, 'c': 9, 'g': 6, 'a': 4, 't': 2, '-:ACTC': 1, '-:ACTCGG': 1, '-:5': 1, '$': 1, '*': 1})

Why is 't': 2' in there?
flags used: -A -B -Q 1 -R -x

Hi Liam, can you show the exact mpileup for this position? Because I see that you tried to parse this into python and also converted the cases (no upper cases).

This way it's easier to see exactly was called.

Thanks @StupidWolf . I stupidly didn't keep the position of the above and it was just one example of a big list. So I'll edit above and give more examples and more details. I appreciate the help. At this point I'm thinking mpileup is a bit buggy, hoping that its resolvable.

Why is 't': 2' in there?

How else do you think the pileup should show you that two reads (in the reverse direction) have a T in that position?

Hey.. the reference is T, so if there are two bases with T on the reverse strand, it should be 2 ","

Yeah, that was the main point.


text:
Are there any tools for discovering tandem repeats around a variant in a VCF file?
I know for example PolyX from Pierre's Jvarkit to find the number of repeated REF bases (homopolymers) around a position in a VCF file as described here exists.
I would like to know however if there tools that describes tandem repeats around a reported variant (TR) up to 6 repeated bases. I can imagine 3 different scenarios:
Scenario 1:
chrK    posY    C     G

TR of (ATC)
____________v
ATCATCATCATCGATGATAC

Scenario 2:
chrK    posY    A     T

             TR of (AACT)
            v____________
CAGTCGATTGGATAACTAACTAACT

Scenario 3:
chrK    posY    C     A

TR of (AG)    TR of (AG)
____________v__________
AGAGAGAGAGAGCAGAGAGAGAG

If you know any tools in mind that can be used for VCFs, I'd really appreciate citing them.
Thanks in advance.

After researching a bit, I found a couple of tools that take as input BAM files. I did not test them yet:

GangSTR
lobSTR
hipSTR
Tandem Repeat Finder (thanks to Michael G.'s answer)

If I find more tools, I will update this post. After testing them, I will come back to share my feedback.

The old program for doing this last updated in 2016 is Tandem Repeat Finder (TRF) here. This will not find a microsatellite sequence at a given variant, but instead find all tandem repeats in the genome and their position. It has been around for a long time.
You would then simply screen this output against the position of your variant. TRF may have been modified however to perform a search within given nucleotide boundaries.

Yes I also saw that software but I didn't give it a proper test. thanks for your Answer !


text:
K means clustering, would PCA be a better option?
I have the data below. I need to use a clustering method to classify them and into categories of "Heterozygotote, Allele 1, Allele 2 and No Call. The values in RFU1 and RFU2 are used to determine the call variable.
Any values in RFU1 and RFU2 that are negative should automatically be classified as No Call. I want to develop my own validation check to check the PCR program calls.
Data
structure(list(Well = structure(1:96, .Label = c("A01", "A02",
"A03", "A04", "A05", "A06", "A07", "A08", "A09", "A10", "A11",
"A12", "B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08",
"B09", "B10", "B11", "B12", "C01", "C02", "C03", "C04", "C05",
"C06", "C07", "C08", "C09", "C10", "C11", "C12", "D01", "D02",
"D03", "D04", "D05", "D06", "D07", "D08", "D09", "D10", "D11",
"D12", "E01", "E02", "E03", "E04", "E05", "E06", "E07", "E08",
"E09", "E10", "E11", "E12", "F01", "F02", "F03", "F04", "F05",
"F06", "F07", "F08", "F09", "F10", "F11", "F12", "G01", "G02",
"G03", "G04", "G05", "G06", "G07", "G08", "G09", "G10", "G11",
"G12", "H01", "H02", "H03", "H04", "H05", "H06", "H07", "H08",
"H09", "H10", "H11", "H12"), class = "factor"), Sample = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
    Call = structure(c(4L, 4L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L,
    3L, 3L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 4L,
    4L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 1L, 1L,
    3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L,
    1L, 2L, 2L, 1L, 1L, 3L, 3L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
    2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L), .Label = c("Allele 1",
    "Allele 2", "Heterozygote", "No Call"), class = "factor"),
    Type = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "Auto", class = "factor"),
    RFU1 = c(-0.295502405, 0.964070798, 3381.332182, 3532.769062,
    3431.836843, 3242.966511, 2104.791167, 2220.008503, 3548.252161,
    3506.51418, 2290.273178, 2281.587684, -5.64819475, -11.73109864,
    3784.914039, 3619.00781, 3618.211608, 3248.106466, 3394.650325,
    3339.870196, 2449.202902, 2426.835174, 3432.153478, 2761.392304,
    -9.267907504, -7.365704356, 3743.092314, 3787.241702, 2172.027787,
    2096.845649, 2135.649551, 2149.145547, 2293.757257, 2348.099108,
    2321.019045, 2022.168867, -17.93532331, -12.59832941, 3805.416768,
    3498.998138, 2304.597239, 2509.63987, 2181.11547, 2261.011876,
    3432.453036, 3662.758933, 2371.11049, 3068.827061, 2612.107589,
    2687.824075, 3179.315918, 3688.525218, 3465.327523, 3405.154043,
    2535.514915, 2452.200079, 374.435299, 423.6015308, 3742.515563,
    3578.777925, 2634.955017, 2527.514043, 3817.579252, 3550.999412,
    -10.72035816, 3294.486334, 3352.40368, 3463.150507, 3472.576514,
    3741.898759, 3571.369947, 3720.645869, 3739.569593, 3855.583168,
    418.6837047, 49.47548241, 2171.034284, 2155.314713, 3432.363384,
    3582.508917, 3425.415274, 3487.203299, 3505.23909, 3413.342096,
    113.5100691, 128.6414053, 2454.588175, 2323.061591, 3188.705702,
    3376.950546, 3291.072437, 3181.001961, 3195.013863, 3776.919197,
    2284.22659, 2277.338631), RFU2 = c(-8.346468029, 235.4058561,
    637.9218251, 650.3759507, 617.4161748, 604.0792911, 4270.310727,
    4199.615749, 689.863543, 712.6144338, 4274.287194, 4541.168491,
    -1.626221758, -2.437395631, 802.0941252, 730.5998997, 686.9037384,
    625.8245403, 644.3586836, 642.8833044, 4937.691887, 5159.479928,
    725.4449756, 573.3910899, -4.006398006, 213.2859144, 739.7910786,
    731.0150586, 4408.81923, 4767.533969, 4302.641493, 4325.913445,
    4597.47663, 4666.904418, 4800.357526, 4142.535329, -17.23239968,
    178.5311942, 778.305843, 743.1438168, 4214.507094, 4553.703511,
    4629.339014, 4459.697405, 661.7299014, 727.1054982, 4553.170272,
    5482.231486, 4520.517999, 4737.802036, 643.3599887, 726.4314715,
    696.5968338, 697.6099599, 411.8118071, 409.4943424, 5687.32635,
    5757.51512, 766.4240193, 779.2403225, 4745.055632, 4582.267792,
    749.5679421, 675.8747055, -7.254521898, 628.3467565, 631.116767,
    672.7064514, 687.2642132, 718.1192917, 731.785499, 668.3686048,
    784.8055727, 791.3155894, 4471.047168, 4501.597841, 4504.670332,
    4442.621066, 682.0632225, 706.6204595, 680.5242182, 683.9558692,
    684.2909706, 618.6535251, 5727.684954, 6098.485474, 5099.952926,
    4779.742057, 571.4303822, 614.9258218, 602.9830491, 651.2847695,
    591.8833499, 742.2387568, 4443.376841, 4716.792177)), class = "data.frame", row.names = c(NA,
-96L))

What I have tried so far
library(cluster)
library(factoextra)
library(formattable)

df <- df[,c(1,5,6)]
df$RFU1[df$RFU1 < 0] <- 0
df$RFU2[df$RFU2 < 0] <- 0
df$RFU1 <- formattable(df$RFU1, digits = 2, format = "f")
df$RFU2 <- formattable(df$RFU2, digits = 2, format = "f")
df$Well <- as.numeric(df$Well)

clusters <- kmeans(df, centers = 4)
Kmeans_plot <- fviz_cluster(clusters, data = df)

This is the plot generated

Points in the top right 57,58,75,76,85,86 should be in another cluster (classed as allele 2 but placed in cluster 1 (rest of points in this cluster are heterozygote)
Also points 24, 55, 56 should be clustered with cluster 2 as all these points are classed as allele 1
I need an algorithm that can detect optimal cluster center, some PCR runs may only have 3 different call variables instead of 4
Programs output - Software company stated their algorithm isn't sophisticated

I would stay away from using k-means, and instead use a method that doesn't a priori define a number of clusters to detect. It also looks as though your clusters aren't exactly spherical, which is an assumption of k-means. I personally am a fan of dbscan, which is available in the R package of the same name. The other poster recommended t-SNE (available in the Rtsne package) but wasn't exactly clear on how / why to use it. t-SNE is a nonlinear dimension reduction technique used specifically for visualization, so you could cluster your points in n-dimensional PCA space and then plot the results in two-dimensional t-SNE space. I think the other poster was recommending plotting your data in t-SNE space before clustering to get an idea of how many true clusters are available, which is another option you could pursue.

Yes, you are right.I will modify the answer

What happens when
clusters <- kmeans(df, centers = 3)

Anyway K means is pretty good and hasn't done a bad job of your data, although I don't know what you mean by the "cluster centre". There are methods to establish 'centroids'.
With K-means you need to vary the group size or have a very good a priori why it is 4, rather than any other number.
The analysis to perform is tSNE ... either by itself or else following PCA. This will provide a clear a priori of the number of groups in your analysis. I don't know R, but it is present in scikit-learn in Python.

Your observation that some of points are in the wrong cluster ... increase the number of groups in k-means, or switch to PCA, tSNE or PCA-tSNE. If the points remain in the "wrong groups" in all analyses at different group sizes for k-means, then basically thats the result. K-means is good, but not perfect, and misclassification will happen.
Another route is the ML repertoire of lasso regression, ridge regression, decision trees, random forests etc ... they all do the same basic job but I'd look at K-means, PCA, tSNE first and if ok ... last. Regression style calculations will give different answers.

t-SNE is a dimensionality reduction approach to produce a bivariate plot and the fondness of this approach is purely observational rather than empirical. It works and also works very well following PCA. Its the method used in Seurat. K-means is used in ML, and certainly has a fan base, but it is not to say it will give you a correct answer.

T-SNE isn't a clustering algorithm, it's a way of taking lots of PCs from PCA and squishing them down so you can see them all in 2 dimensions.  It can help you eyeball clusters which will be based on more than just the top 2 PCs.  But it won't literally identify clusters or assign points to clusters.

Oh of course I know clustering statistics very well because they were originally used in phylogeny. The problem with clustering algorithms is their is a lot of them - they give different answers - but k-means has gained respectabilty within ML. However the OP wants "clusters" in a non-statistical sense - groups/discrimination/centroids

I don't think clustering algorithms are working here, it looks like you can do this much simpler by looking at whatever determines the X axis in your plot.  The samples split cleanly when looking in just one dimension.

Hi @swbarnes2 thank you for your opinion. Could you provide the code you used? I have been trying to work this out for a while.

Why are you making well numeric?  Isn't that going to make it be included in the clustering calculations?

the well variable was originally a factor, and k-means only runs on numeric. I just want to develop an algorithm that can determine how many centroids per dataset. I have multiple PCR outputs and have to run this algorithm to verify if the calls generate from the program's algorithm are accurate. (I have attached a screenshot of what the program generates. If I plot the RFU1 and RFU2 how can I use clustering to classify and return "no call, allele 1, allele 2 or heterozygote"

Then you should remove the well from the data columns, not make it a number that gets incorporated into the kmeans clustering.  You have 2D data.  You can't run PCA, let alone t-SNE on 2D data.  You don't need sophisticated algorithms to handle 2D data.

I used the elbow method to find the optimal cluster numbers and ran k-means following this. Every time I run the k-means I get different clusters. How can I run that I get the same cluster pattern and cluster numbers each time?

Honestly, I would stop trying to do this with pure math.  Just use your RFU values, and your eyeballs, and divvy up the samples by eye.


text:
Drawing 2 heatmaps with the same gene order
I have 2 dataframe with same genes both different data sets
like
> head(heat1[1:3,1:4])
             A10         A11        A12        A2
ATP5F1 2.1163092  0.32044094 -0.5661269 0.7978236
DDX5   0.6277182 -1.67055853 -0.0123082 0.8895472
EEF1G  1.6459212  0.02923401 -1.2737676 0.5842162
>

> head(heat2[1:3,1:4])
              A10        A12          A3          A4
A2M    -0.6031969  0.6378864 0.000573381 -0.07769313
ABCB11  0.5600164  0.3087120 0.171636836  1.25681512
ABCC2   0.2093133 -0.5608016 0.782278709  0.16618683
>

I want to plot them side by side with the same order of genes
I tried
> ht1 = Heatmap(heat1)
> ht2 = Heatmap(heat2)
> ht_list = ht1 + ht2
> draw(ht_list, ht_gap = unit(1, "cm"))
Warning message:
Row names of heatmap 2 is not consistent as the main
heatmap (1)

But I getting not the same order of genes

Can you help me with this?

Please read through the ComplexHeatmap manual, especially the section on row and column orders: https://jokergoo.github.io/ComplexHeatmap-reference/book/a-single-heatmap.html#get-orders-and-dendrograms-from-heatmap
You can use row_order and column_order following a draw() (which you can nest between a pdf(null) and dev.off() to disable drawing it on any device). That way, you'll be sure to get the right order from the first heatmap to use with the second.
By nesting, I mean like so:
h1 <- Heatmap(...)
pdf(NULL)
h1d <- draw(h1)
dev.off()
ro1 <- row_order(h1d)
co1 <- column_order(h1d)

Use whichever (row/column) order is appropriate to gene order.


text:
Tools for comparing/visualizing FASTAs?
I have two FASTAs/assemblies of the same species, but the bases are somewhat different. I would like to explore this.
What tools/methods exist to compare two FASTAs and see the difference in structure between them? Are there visualization tools available, something like IGV?
I think users would need to do a multiple-sequence alignment between these FASTAs, and then check for the differences between them. However, I don't know if a "standard approach" exists to do this.

You can try Mauve multiple genome aligner.

Mauve is a very good tool for this! Thanks. Does Mauve allow for us to visualize base-level differences (easily)?

I am not sure. I have not tried it at that resolution and i have not used it in a long time.

Depending on how long the sequences are, you could use the LAST Web Service:
http://lastweb.cbrc.jp/
It has an HTML view that shows base-level differences with colour-coding.
LAST [download / command-line version] can also be used for producing whole-genome comparison dot-plots.
Another thing I do for comparing assemblies is to chop one assembly into smaller pieces (e.g. 10kb), map the pieces to other other assembly to generate BAM output, then look at the BAM pileup file in a BAM viewer (e.g. Tablet).

For whole-genome align:
Maybe you could use mummer? Mummer is good to align the whole genome of species, and it also has so many other tools to modify the result. (It can show the SNP.)
For short sequence:
You could try mega X

Before you compare your genome assembly fasta, I would recommend you first do a reference guided scaffolding using the tool RagTag.
Once you are done you can compare the newly generated scaffolds (with other genomes or not) using the software BRIG. I made a video of how to use this software here
In addition you can also do a genome annotation on your scaffold using prokka. I made a video for that as well here
Good luck


text:
How are these definitions related to differential expression?
I have two groups of patients; for each patient I have an output file (RNA-seq) contains this information
                  readcounts_union  readcounts_intersectionNotEmpt  genelength  R/FPKM_union    R/FPKM_intersectionNotEmpt
ENSG00000000003.10             8                      8               3688     1.06120392730222        1.05722142116966

ENSG00000000419.8             31                     31               4032     3.76132572546529               3.74721015127372

People say
readcounts_union

reads that overlap a single feature are included in the count for that feature; reads that overlap more than one feature are not included in the counts
readcounts_intersectionNotEmpty

regions that are shared between features are excluded; reads that overlap a remaining region are included in the count for the feature that includes that region; reads that overlap regions from more than one feature are not included in the counts
R/FPKM_union

readcounts_union normalised for gene length and sequencing depth
R/FPKM_intersectionNotEmpty

readcounts_intersectionNotEmpty normalised for gene length and sequencing depth
I want to do differential expression analysis between groups but I don't know which column I should use for DESeq2
Anybody familiar with such a format?

I've never seen this format before, what program generated this file? DESeq2 requires raw counts as input (it does it's own normalization internally), so you should avoid using the FPKM values with DESeq2

Your life would be easier if you simply used standard formats instead of whatever this is. Literature currently seems to prefer pseudoalignment and/or selective alignment (kallisto, salmon) as quantification method, followed by gene level aggregation and diff. analysis. Why bothering with such non-standard things? Get the raw fastq files and start over. salmon-tximport-DESeq2, fully covered with example code in the DESeq2 vignette and tutorial at Bioc.

I don't have access to the raw files, I rather have been given only this format and been asked to give them list of differentially expressed genes

Either of the first two columns are fine to use, though I personally prefer the "intersectionNotEmpty" one.
As an aside, that looks like a custom format file produced from post-processing htseq-count output.


text:
What samples can be used for a Mutect2 Panel of Normals
I'm working on calling somatic variants from solid tissue tumors. I have one normal and one tumor sample for each patient. I plan to use Mutect2 to call somatic variants after preprocessing the data with the GATK Best Practices pipeline. Mutect2 recommends the use of a Panel of Normals (PON), but I can't figure out what samples should be used to built the PON. The guide says

Additionally, the samples should come from subjects that were young and healthy to minimize the chance of using as normal a sample from someone who has an undiagnosed tumor. Normals are typically derived from blood samples.

My normal samples are from solid tissue, not blood.
Can I use the normal samples from each patient to build the panel of normals? Or would this cause problems because I'm using reusing the same normals for both the PON and as the normal sample for somatic variant calling? My first thought is that most somatic variants in individual normal samples should be filtered out by pooling all the normals together to build the PON but I wasn't sure if there are other reasons to avoid using non-blood based normal samples to build the PON

It just says normals are typically derived from blood samples. It does not say there is something wrong with not using blood.

As suggested by @burger in their comment, you can indeed use solid tissue samples as normals for your PON. You are clearly aware, but I want to spell this out for others, that there are a few caveats that go with the practice of using solid tissue normals rather than blood normals. Mainly, if your normals are "tissue-adjacent", ie sampled from an area of tissue adjacent to the tumor, there is a higher risk of including contaminating tumor cells compared to blood normals. As a result you might have real somatic variants that end up filtered because they are seen in the normal. Though that typically acts more at the level of the direct Tumor-Normal comparison than at the level of the PON, since the PON draws from many normal samples and the chance of seeing the same contaminating somatic mutation in many such samples should be quite low. So ultimately that should be fine.
It's a very common practice to use whatever matched normals you have available to build your PON, and there doesn't seem to be any drawback to doing that as long as you have enough samples (at least a few dozen, ideally a few hundred or more). If you don't have enough samples in your own study, then you should seek to use other samples from young and healthy people that have been sequenced using the same technology (same library prep kits/protocols, same chemistry version etc) since the goal of the PON is to weed out recurrent technical artifacts.


text:
gnu parallel for macs peak calling
These are the list of files
Control_Input_sorted.bam_rem.bam
Control_H2BUb_sorted.bam_rem.bam
Control_IgG_sorted.bam_rem.bam
PTPN6_g2_6_H2Bub_sorted.bam_rem.bam
PTPN6_g2_6_Input_sorted.bam_rem.bam

So to get unique files i do this
 ls -1 *bam | sort | sed -r 's/_sorted.bam_rem.bam//g' | sort | uniq
Control_Input
Control_H2BUb
Control_IgG
PTPN6_g2_6_H2Bub
PTPN6_g2_6_Input

Now i have to run each sample against input and IgG for peak-calling.
Like  Control_H2BUb_sorted.bam_rem.bam against
Contol_Input_sorted.bam_rem.bam & Control_IgG_sorted.bam_rem.bam and

PTPN6_g2_6_H2Bub_sorted.bam_rem.bam against PTPN6_g2_6_Input_sorted.bam_rem.bam

So how do I parse and execute using gnu parallel.
I came across this tutorial where it does something like this
cat sample_names.txt | parallel --max-procs=12 'macs2 callpeak -t {}-A-NC.sorted.bam \
 -c {}-G-NC.sorted.bam -g hs -n {}-A-NC-sharp-model -q 0.01 --outdir {}-A-NC-sharp-model-peaks 2> {}-A-NC-sharp-model.stderr'

Now im not sure how to pass the file names as argument for IP and Input or IgG.
If i use logic a little bit i have to fix Control and PTPN6_g2_6 as constant naming!.
Any suggestion or help would be really appreciated.

not gnuparallel but nextflow. I don't know macs and what is 'first.bam' and 'second.bam' is not clear not me but I hope you'll get the idea.
The nextflow:
Channel.fromPath(params.bams).splitCsv(header: false,sep:'\t',strip:true).map{T->file(T[0])}.into{bams1;bams2}

boolean isControl(file) {
    return file.name.endsWith("_Input_sorted.bam_rem.bam") || file.name.endsWith("_IgG_sorted.bam_rem.bam ");
    }

bams1.filter{T->isControl(T)}.set{controls}
bams2.filter{T->!isControl(T)}.set{experiment}

process mac {
    maxForks 14
    tag "${bam1} vs ${bam2}"
    input:
        set bam1 , bam2  from experiment.combine(controls)
    output:
        file("output.peaks") into mac_peak
    script:
    """
    # check syntax please
    macs2 callpeak -t {bam1.toRealPath()}  -c ${bam2.toRealPath()} -g hs -n out-model -q 0.01 --outdir output.peaks
    """
    }

invoke:
find /dir1/dir2 -type f -name "*.bam" > input.list

nextflow run --bams input.list workflow.nf

okay let me give it a try first and update you.

my first try ended up like this "Not a valid project name: workflow.nf" I need to read the basics of nextflow

first bam is my Control_H2BUb and PTPN6_g2_6_H2Bub , and it should be compared with samples labelled Control_Input,Control_IgG and PTPN6_g2_6_Input which are second bam. But the catch is the Control lablled samples should always be compared with the control labelled and same goes with PTPN6 labelled


text:
Sequence alignment of extended IUPAC sequences
I want to align a short sequence to a reference, and both are encoded in the extended IUPAC notation (i.e. 16 letters).
What alignment software can I use for that?
For example, in Bowtie "all non-A/C/G/T characters in reads become Ns", so it ignores which extended letters are matching.

Found a very similar question here: https://biology.stackexchange.com/questions/7919/tool-for-nucleotide-alignment-with-all-nucleotide-codes-e-g-r-y-w-s-etc

You can try using Biostrings, depends on how large is your reference, if it is small, you can just pairwiseAlignment with the nucleotideSubstitutionMatrix which supposedly :
library(Biostrings)
data(phiX174Phage)

x = DNAString("TCCAWRSCGCAGAAGT")

lapply(phiX174Phage,function(y)pairwiseAlignment(x,y,
substitutionMatrix = nucleotideSubstitutionMatrix()))

$Genbank
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$RF70s
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$SS78
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$Bull
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$G97
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$NEB03
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

If it is much larger, you can convert the sequence of interest to a pwm and find it in the reference fasta, like shown in this example


text:
Economist article on coronavirus
I am wondering about an article in the Economist here:
https://www.economist.com/briefing/2020/02/29/covid-19-is-now-in-50-countries-and-things-will-get-worse
There is a graph there

The explanation is as follows
The course of an epidemic is shaped by a variable called the reproductive rate, or R. It represents, in effect, the number of further cases each new case will give rise to. If R is high, the number of newly infected people climbs quickly to a peak before, for want of new people to infect, starting to fall back again (see chart 2). If R is low the curve rises and falls more slowly, never reaching the same heights. With sars-cov-2 now spread around the world, the aim of public-health policy, whether at the city, national or global scale, is to flatten the curve, spreading the infections out over time.
I am having trouble with this. There are 2 curves here. First of all, I am not sure of the different processes that determine them. It says one is without measures. The other is with measures. But it also says one is with high R and the other with low R. Is R a fixed constant for a disease? Alternatively, is there some variability in R and does that degree of variation affect transmission? This is not explained or stated in the article and I haven't heard R described in this way. At best it is not written clearly in my opinion.
In addition, does without measures necessarily imply high R? Is that what they are getting at? Also, it looks to me as though the blue graph goes to zero very quickly, basically when t is equal to the mean of the orange curve. That means that without measures the epidemic will end quickly (because it exhausts all the possible patients).
Finally the areas under the two curves is about the same. This means that the number of patients is about the same with or without measures.
Anyway, can anyone clear up any of this?

R is a reference to R0, otherwise known as the basic reproductive rate, and means the number of new cases from a single patient infection. R0 has to be above 1 for a disease to persist. It is a classic equation which account for the probability of transmission against immunity. The y-axis is the number of infecteds, very poorly presented.
The idea is that R0 will be lower if there is containment (less people to infect because you can't move around), but will last longer, possibly because the number of naive patients (not encountered the virus) will be larger over a longer period of time.
That particular graph appears to be saying that everyone at risk of infection will get infected with or without control measures (the surface area under the graphs is approximately the same*) and the only thing that will halt the virus is acquired immunity. If this is mainstream economic thinking then I beg to differ, albeit I note the source is CDC (again see note*).
For example, all other respiratory viral infections are seasonal. If COVID-19  is seasonal containment is the right choice, because the window of transmission is limited.
*, Given this is of CDC origin I think either its a mistake, or at the very least was taken out of context, because the idea is that overall there is a lower number of people infected as a result of "social distancing" and that does not appear to be what the two curves show.

For example, all other respiratory viral infections are seasonal. If COVID-19 is seasonal, containment is the right choice, because the window of transmission is limited. This part is valid.

That particular graph appears to be saying that everyone at risk of infection will get infected with or without control measures (the surface area under the graphs is approximately the same*) and the only thing that will halt the virus is acquired immunity. THIS Part is not in consonance  with the main theme and interpretation of presumptuous normal  curve with greater height can not be substantiated.

Hi @subhush, whilst I agree the curve needs integrating and is at present a cartoon, it doesn't look optimistic.

R, the reproductive number, relates to the average number of (new) people that will get sick (infected) per person that is already sick.
For instance, if $R=2$ and you start with a single infected person, then the next generation will be 2 people, those 2 people will make 4 people sick, those 4 people will make 8 people sick, and so on.

Is R a fixed constant for a disease?

This reproduction number, how many people the disease is passed on to, depends a lot on the circumstances. For instance, how frequent are contacts and how likely is transmission per contact. The probability of transmission per contact will depend on the type of disease (and from that stems those comparisons of different disease), but can also be influenced by measures like reducing frequency and intensity of contacts (social distancing).
The coronavirus spreads in droplets from breath, coughs and sneezes, and will reach other people mostly during close contact (Although there may be events of long-distance spreading if the virus is spread in aerosols, which is claimed to be not the case but it may still exist to a small extend, like with SARS in 2003 explained here and here). So, the reproduction number, R, will likely be influenced by measures like preventing close contact, what they also call social distancing, as well as proper hygiene (washing hands to prevent transmission from virus picked up from fomites).
A clear example is the case of a sexually transmittable disease, where condom use, or not having sex at all, will reduce the spreading of the disease.

Alternatively, is there some variability in R and does that degree of variation affect transmission?

There is a lot of variability in transmission rate.
The spread of the covid-19 virus has likely been accelerated by celebration of New Year's in China and skiing holidays and carnaval in Europe (as well as large scale meetings during soccer matches and church services).  (this is why large events are bad)
Also, the disease will at some point satiate the surrounding. Once people have been infected they are likely immune for some while and less likely to be infected again. So once the entire household, classroom, or workplace is infected then new infections will be less likely. (this is why long distance travelling is bad because it brings the virus quickly into untouched area's)
It is wrong to speak about the reproduction number. The reproduction occurs heterogeneously and it is wrong to speak of a single number.
The average reproduction might be 1 or less, but there can still be effectively exponential growth within a sub-population with relatively higher transmission rates.
Especially when this sub-population is strongly connected network on a large scale. (the virus is like hopping from cluster to cluster in a large network, the transmission occurs mostly in dense clusters of the network and some area's remain relatively little affected, this is why Sweden does not need to take though measures, their network is relatively sparse and safe)

Finally the areas under the two curves is about the same. This means that the number of patients is about the same with or without measures.

That sketch in the economist article is just a simplistic sketch. The 'flattening the curve' will reduce the rate of growth, but it should also reduce the total number of cases.
An example is the image below from this SE post:
https://stats.stackexchange.com/questions/461319/using-different-nodes-in-a-networked-compartmental-model-sir-for-different-reg
In the image you see the effect of different social distancing measures (more details about it in that post). On the left the cumulative cases are plotted (this relates to the area under the curves) and you see a reduction of the total cases. The curves are not only flattened, but also reduced in volume. (Of course, this does require that the measures are continuously followed for the entire period. That image is an idealistic situation)

Also, it looks to me as though the blue graph goes to zero very quickly

The endemic will eventually halt when a lot of people have already gotten sick (it may also continue for a long time, and become endemic, because the population is not static, and births introduce new people that have not been sick before, also immunity may become lost over time).
The basic reproduction value, often denoted $R_0$, is the number of people that the infection is passed on to when nobody in the surrounding is immune. But if a lot of people have gotten sick already, or when there is a vaccine, then a lot of people are immune, and less new people will acquire the disease. Effectively the reproduction will be less.
When the disease spreads faster, then this point of satiation/immunity will also be reached more quickly. Currently, strategies are being developed to make the disease spread as fast as possible (to get over with it quickly) but not too fast (to not overload the intensive care units and health care system). This may take a long time (there is still a lot of uncertainty about it) and possibly not economically viable.
If you see that the spreading might be over more quickly when we do nothing about it, then it could be better to just let the disease continue without measures, accept the high amount of deaths, and be done with it quickly instead of letting it linger for a long time. (But not many would dare to suggest these things seriously, because of the complicated ethical issues, and a stigmatization of an utilitarian approach to health and preventing deaths. )


text:
Scooping and GeneMatcher
A colleague has uploaded about a dozen candidate disease genes to GeneMatcher, described here and here and located here,  with excellent success. However last year and again recently he has encountered the case wherein a case was uploaded and a cohort formed, but a few months down the line one of the dropped out collaborators published an n=1 paper in a low tier journal, basically using GeneMatcher to check whether the disease link of the gene was legitimate, to assess the phenotypic range and see what test are worth doing.

Is there documented precedence of this behaviour (preferably citable)? If not,
does anyone have personal experience of this and what is the ratio of incidences over total submissions (I know that is hearsay, but it would be very helpful)?
And also has anyone successfully taken action (say getting user barred)?

when you say "this behavior", do you mean an ex-collaborator trying to rush a part of the collaborative project into print to get credit? that is very common. Are you referring more narrowly to this behavior on the GeneMatcher platform?

I was not aware that is common in general for a collaboration to be formed only for the collaborator to drop out and rush to press with the intel. Turns out there are lots of collaboration agreements to protect against shenanigans in joint works which are similar to MTAs. Sickening as it sounds this is the solution.

I think that often collaborator starts with good intentions, but then the team falls apart and stop liking each other, and they still want to get something out of the work. There are a few ways it can happen. Never happened to me personally (at least not like this specifically). But anecdotally, I'd guess this is actually one of the most common paths to scooping.

Ah, no, this is more malicious. The site matches geneticists and they exchange patient details in chain emails, discuss a path to publication (agreed by all) and make the clinicians request further analyses say MRIs or blood work (within 2-3 months). Then when follow-up lab tests are planned, the collaborator suddenly goes silent and half a year later there is a paper in a low tier journal with their sole patient but with strong emphasis on the key traits discussed on in the collaborations.

Ah, I see. that is a more specific and very bad behavior.

I am having trouble finding anything that directly addresses the issue of failed collaborations leading to scooping, but here are some examples of similar behavior:

it is sort of implicit in the discussion here.

This question on academia SE is related but not exactly right.

There are some related stories here.

Here is an article about plagiarism among collaborators, which is quite similar.

I feel like a lot of the discourse around preprinting etc. is similar but again not quite there.

As for punishment, that sounds hard. Scooping is sort of built into science; fighting over precedence is a big deal and has been for a long time. My solution would be to get rid of current academic publication and have everyone's open digital lab notebooks streaming results to publication feeds in real time; that way there's no ambiguity. Preprints and micropublications are sort of a poor man's version of this.
This is not a very complete answer but there may be enough there to help.

Thanks, excellent points. Just to add some extra in terms of GeneMatcher. Most collaborations flowered and only a few sour. Someone else has had a similar experience to my colleague, so it is not uncommon. The solution is to use collaboration agreements to safeguard against this —even if these are laughably enforceable.


text:
Using cellranger mkgtf w/ NCBI GFF3 files
I am using several cellranger tools/pipelines to analyze some data, specifically, mkgtf, mkref, count, and aggr. I am using mkgtf and mkref to include a PI-synthesized protein to my reference. I have finished the analysis but am not seeing the protein in my output, so I am going back to the start to verify the process.
Unlike the documentation (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references), my GTF file was not downloaded from ENSEMBL or UCSC, but NCBI (https://www.ncbi.nlm.nih.gov/nuccore/AY678269.1?report=genbank). This GFF file does not have the same key nomenclature as seen in ENSEMBL or USCS.
I did use the mkgtf tool, but the created file is identical to the input.
cellranger mkgtf custom_ref/AY678269_genes.gtf reference_data/AY678269_mkgtf_genes.gtf

I went ahead and used cat to add the fasta and gtf file from the synthesized protein to the respective reference fasta and gene files, then ran mkref.
The s.protein is in the fasta file:
>lcl|AY678269.1_cds_AAV52169.1_1 [protein=tandem-dimer red fluorescent protein] [protein_id=AAV52169.1] [location=1..1431] [gbkey=CDS]

as well as the gene list
##sequence-region AY678269.1 1 1431
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=32630
AY678269.1      Genbank region  1       1431    .       +       .
AY678269.1      Genbank CDS     1       1431    .       +       0

However this protein is not found in the output...surely I must have made an error. Does anyone have any ideas about what I could be doing incorrectly?
Do I need to change the entire nomenclature of the GFF keys to match what I see from ENSEMBL (e.g., gene_id ""; gene_version ""; transcript_id, etc...)?
Any assistance is greatly appreciated.

Regarding your comment:

This GFF file does not have the same key nomenclature as seen in
  ENSEMBL or USCS

The link you have provided corresponds to an entry in Genbank format and not to GFF format. I believe you need a file in GTF format in particular.
Nevertheless, given the fact that you have manually added the corresponding entries to the reference fasta and gtf files prior to mkref, you should have seen your gene of interest in the output. I would double check one of the original fasta headers and original gtf entries in the reference file and match these for the new entry.

Thank you for the response. I figured out what the problem was, I was given an incorrect accession number! However, after finding & reading through this: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/tutorial_mr, I can see how my approach may have been problematic anyway because of the NCBI keys. Fixed it & I'm creating the custom reference. FYI, the link does go to a genbank file, but when downloading the information, you have the option of choosing the format (fasta, genbank, gff), but you are absolutely correct about the GTF vs GFF. Danke!


text:
Filtering raw sequencing reads
I have to filter the raw sequencing reads based on the following criteria:

Remove reads containing adapters
Remove reads containing N > 10% (N represents base that could not be determined)
Remove reads where the Qscore (Quality value) of over 50% bases of the read is <= 5

The reads are paired-end fastq files. I have been able to do the first part by using trimmomatic:
trimmomatic \
  PE in_1.fq.gz in_2.fq.gz \
  pair_out_1 unpair_out_1 pair_out_2 unpair_out_2 \
  ILLUMINACLIP:adapter.fa:2:30:5

The other 2 steps are a bit more challenging since I have not found a tool that perform exactly what is required. I could write some custom code of course but I would like to use the existing tools if possible. No need to reinvent the wheel.
Do you have any suggestions on how to proceed?

I think you can come close to this with cutadapt:
cutadapt \
  -a ${FWD_ADAPTER} -A ${REV_ADAPTER} \         ## adapter
  --max-n 0.1 \                                 ## do not allow > 10% Ns
  -q 5 \                                        ## remove bases with basequal < 5
  --minimum-length (readlength/2) \             ## minimum length after trimming
  -o read1_trimmed.fq -p read2_trimmed.fq \     ## outputs
  read1.fq read2.fq                             ## inputs

This will remove reads with more than 10% of bases being N and remove bases with base quality below 5. If you then set --minimum-length to half the read length, this should also remove reads where more than half of the bases had a base quality below 5 and therefore got trimmed away. Poor base quality in Illumina reads typically happens towards the 3' end so I think this little workaround is a sufficient proxy to achieve your goal. A base quality of 5 is very low in any case, so if sequencing went normal the typical quality should be notably higher. I would be surprised if this third filter is really necessary and will eliminate a notable number of reads.


text:
How to run trimmomatic in HPC
Sorry if this has been asked before, I've done a quick search and I don't think there's an easy explanation for me to understand.
I'm really new to bioinformatics/RNASEQ analysis, having only taken up an attempt to process my RNASEQ data last week. I have 0 background and some rudimentary knowledge of command line. So far I've managed to get access to my universities HPC and perform a fastqc on my fastq data.
What I'm stuck at at the moment is attempting to trim my RNASEQ data and remove adapters. My data had been run through Novaseq and I'm not exactly sure how to trim my current RNASEQ data.
I've attempted to run the following script using trimmomatic
# go to working folder
cd /project/CHR_RNASEQ/reads

# run trimmomatic
trimmomatic SE -threads 8 -phred33 SLIDINGWINDOW:4:20 ILLUMINACLIP:$TRIM/adapters/TruSeq3-SE-2.fa:2:30:10 MINLEN:36 1_R1.fastq.gz 1_R1_trimmed.fastq.gz

It did log a job in the HPC, but I had no output at the end. I didn't receive an error either. This is unfortunately just some code I had lifted and trialed and errored until something finally ran, but I'm very sure I'm very wrong and any help with understanding this would be really appreciated.
Cheers

If you're using a HPC please talk with your IT department to check which job scheduler they are using: LSF, SGE, SGE, SLURM, PBS. Most of cases you are not able to run any job without this queue system.

Depending on the cluster management tool, you might have received e-mails when the "job" begins and ends. If so, you can check the "Exit status" of the job. For example, in the case of our HPC the relevant lines from the "job completion e-mail" would be:
Execution terminated
Exit_status=0

BTW, at first I was confused with the phrase "Execution terminated", however, in the scheduler used in our system, it means "the run has been successfully or unsuccessfully completed". Exit_status=0 int the next line means a successful run and for the unsuccessful runs there is a list with the error codes.
Secondly, if your run was successful, you would see new fastq files containing "_trimmed" in the filename (reading this from your code). These fastqs should be more or less similar in size to the original depending on how aggressive the trimming settings are.
If you can provide more information, i.e. the type of cluster management tool used, the log file, ..., you can get more help and precise answers from the community.

okay, I'm really sorry, but I'm not sure what you mean by cluster management tool or log file. I tried looking around for information on this, and did get kind of stuck. I looked in the HPC for my file, and couldn't find one, so i'm assuming it ran the job but didn't actually give me an output.

You wrote "It did log a job in the HPC", what is in this? In my case, there are two logs at the folder from which I submit my job, one like "jobname.ejobid"  and the other "jobname.ojobid". The first is the error log and the second for the output from the software you are running. Do you have these?