Since there doesn't seem to be an easy way to run this in parallel, you could instead break the job into sections. For example, separate each chromosome into its own gtf file, extract the sequences using that file and then cat them all together.

The commands would be something like this (using the human hg38 genome and the gff GENCODE annotations downloaded from here):

Split the file into one file per chromosome (or whatever else you have as the first field of your gff/gtf file):

tmp=$(mktemp -d); ## create a temp directory
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf

Extract the sequences described in each file, running each command in the background (&):

for f in $tmp/*; do
    gffread -w "$f.fa" -g hg38.fa "$f" &
done

Concatenate them into one file:

cat "$tmp/*fa" > all.fa

You can combine all this into a single command, using wait to make sure you don't concatenate until all of the sub-commands have finished:

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf &&
    for f in "$tmp"/*; do
        gffread-0.9.12.Linux_x86_64/gffread -w $f.fa -g hg38.fa "$f" &
    done
    wait
    cat "$tmp"/*fa > all.fa

Or, just save it as a script called pgffread.sh (or whatever):

#!/bin/sh
gtf="$1"
genome="$2"

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' "$gtf" &&
    for f in "$tmp"/*; do
      gffread -w "$f".fa -g "$genome" "$f" &
    done
wait
cat "$tmp"/*fa
rm -rf "$tmp"

And then run it giving the genome and annotation files as parameters:

pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa

On the system I ran it on, which has 48 cores, the time difference was significant (although it's so fast anyway, I really don't know if it is worth the effort unless you do this very often):

$ time gffread -w out.fa -g hg38.fa gencode.v27.annotation.gtf
real    1m52.072s
user    1m3.160s
sys     0m47.257s

$ time ./pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa
real    0m17.533s
user    1m5.228s
sys     0m43.039s


how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id)

If you can use awk, you could use this simple script:

awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^>/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta

You can also put the code inside a script file:

extract_genes.awk

#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^>/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?

And call it like this:

awk -f extract_genes.awk gene_list cds.fasta


FeaturePlot from Seurat: change its title

I tried with some data that I have and this is working for me:

p <- FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                      reduction.use = "tsne", do.return = TRUE)

lapply(p, function(x){x + labs(title = endothelial_symbols[1])})

I think it is because FeaturePlot returns several ggplot objects according to the ids you put. This is returned as a list of plots so you have to iterate each one for labeling.


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.


Since Spearman is a rank-based test, it relies on you being able to accurately decide on the ranking of your observations by some metric (usually the magnitude of the numbers). If two observations have identical values (are tied), then they cannot be definitively ranked. Since the ranks are not unique, exact p-values cannot be determined.

e.g. in your data, for gene 2:

                  gene2 rank
TCGA-2Y-A9GS-01A    7.4    1
TCGA-2Y-A9GT-01A    1.0    5
TCGA-2Y-A9GU-01A    2.0   3=
TCGA-2Y-A9GV-01A    0.5    6
TCGA-2Y-A9GW-01A    2.1    2
TCGA-2Y-A9GX-01A    2.0   3=

By way of example, it's possible to do the test with either TCGA-2Y-A9GU-01A or TCGA-2Y-A9GX-01A ranked number 3 - meaning a definitive p-value cannot be determined.

If there are only a few ties, you probably don't need to worry too much. If your FPKMs (or better, TPMs) were calculated with greater precision, then ties would be less likely. Is there a reason you only have these numbers to 1 decimal place?

Log or linear scale doesn't make a huge difference (and neither is 'wrong'). If you choose to plot log counts, you need to consider what happens if you have zeroes in your data.


Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis

It is indeed true that for the DE analysis one should include batch into the formula to avoid changing the original counts. Still, for everything else such as plotting heatmaps use of removeBatchEffects is perfectly fine and (at least to me) a standard and well-accepted procedure. It essentially does not matter what you use to correct for the batch effect for the counts you use downstream. Results will probably be similar. Combat-Seq from the sva package is a recent adaptation of ComBat specifically for RNA-seq which (from what I understand) better deals with the integer-count nature of the data. This operates on raw counts and avoids the infamous negative values that happen at times with both limma and Combat. After applying ComBat-Seq on your raw counts you can normalize them as usual with edgeR (or any tool you like) and then make the heatmaps. See https://github.com/zhangyuqing/ComBat-seq. The DE results still should come from the normal DE pipeline with batch as covariate as discussed above.


Merging files gives memory allocation error

You could try something like this :

First you should work with the same class of object : data.table for fast merging operation

# Turning a to a data.table
setDT(a)
# b is already a data.table

Then you might have to set the keys of both table :

# I added chr in case you need to match the same position on a different chromosome.
setkey(a, chr, position)
setkey(b, chr, position)

Finally perform the joining operation :

a[b, nomatch=0]

Because you didn't provide the expected output, it might need some tuning !


STAR quantMode (GeneCounts) essentially provides the same output as HTSeq-Count would, ie. number of reads that cover a given gene. This is the most simple measure of expression you could get from RNA-seq data. Kallisto and Salmon utilize pseudo-alignment to determine expression measures of transcripts (as opposed to genes). RSEM uses some algorithm to determine isoform fractions, but still based on alignments. (Note that I don't have any personal experience using RSEM.)

It ultimately depends on the research question which method will be most appropriate for you.

read counts

pros:

relatively easy to understand concept-wise
can be used to generate many different expression measures (FPKMs, TPM, CPM, etc.) and is, thus, more versatile in down-stream analyses

cons:

you can't distinguish isoforms (from the counts themselves)
requires alignment and, thus, more time to obtain

salmon/kallisto:

pros:

faster (due to lack of alignment step)
quantifies (known) isoforms

cons:

expression values are more abstract, less easy to explain
gene-level analyses will require some extra work

RSEM

Likely lies somewhere in between, still requiring alignments will make it slower. But it does give you expression values for isoforms.


In stageR, we automatically return stage-wise adjusted p-values of zero for transcripts coming from genes that were found to be significant and only have 1 transcript, when testing for differential transcript expression (DTE).
Similarly, when testing for differential transcript usage (DTU), we return p-values of zero for transcripts derived from a significant gene consisting of no more than two transcripts (based on the tx2gene object you provide).

Transcript-level p-values for non-significant genes should not occur. We provide NA's for genes with only one transcript for DTU, and by default you should get NA p-values for non-significant genes, since these genes do not proceed further than the screening stage.

Could you provide some more details on the analysis? Which multiple testing correction have you selected?

Would it be possible to provide a reproducible example for your analysis?


See IUPAC codes:

So, as you can see above, K means "Either G or T".


What are the units for the force output by nstfout, the velocity output by nstvout, and the coordinates output by nstxout in gromacs?

nstxout, nstvout and nstfout are the number of steps between coordinates, speed and forces written in the output trajectory (Output control documentation).
As for the units present in your trajectory files (say a .trr):
coordinates are in nm, speeds in nm.μs⁻¹ (so in mm.m⁻¹) and forces in kJ.mol⁻¹.nm⁻¹. These should be consistent through all GROMACS files (Definitions and Units documentation).
Note that .xtc files multiply coordinates by 1000 to compress them, so they are stored in pm (see GROMACS xtc documentation). Given MDAnalysis documentation on the subject, I believe .xtc files can only store coordinates, no speeds or forces.
I would be wary of softwares used to analyze your trajectory results (MDAnalysis, VMD, etc), as they may automatically do the conversion of these units, and for example show you coordinates in Å.


What database do gene IDs starting with sp_v2 refer to?

This is possibly a wild goose chase, but a lot of searching led me to this sample:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2536852

Which appears to exhibit accessions of the right format. This in turn leads to the platform:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL22166

Which is a custom spotted cDNA array, which the metadata seems to suggest is a human platform.

The full annotation table gives gene symbols and Entrez IDs for most of the probes on the array:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL22166&id=53508&db=GeoDb_blob144

To take the example of the first line of your results:

ID          ORF     Entrez gene SEQUENCE
SP_v2 4634  ETNK1   55500       AAAGCAGCTTCATCTTTCAAAATTGATTTGCTCTGGTTTT

The Entrez Gene record for that ID matches up:

https://www.ncbi.nlm.nih.gov/gene/?term=55500


Testing for viral/bacterial infection FASTQ files

If you're talking RNA, I have a blog post from a while back about Diamond/MEGAN:

http://achri.blogspot.com/2017/04/data-for-free-viral-metagenomics-in.html

In principle this will work for DNA if you have a non-trivial number of pathogen reads (since they tend to be gene dense in their genomes). Otherwise I use Kraken on a database for all the NCBI nt database, but this is slower and doesn't let you sanity check the false positive taxonomy assignments as well as I describe in the post. Of course, ideally, you run multiple methods if you have the resources to do so.


Within and between sample count normalization

The approach you are describing seems very strange.
Crucially, the Vignette for DESeq2 states that the model only works correctly with unnormalized counts as input:

It is important to provide count matrices as input for DESeq2’s
statistical model (Love, Huber, and Anders 2014) to hold, as only the
count values allow assessing the measurement precision correctly. The
DESeq2 model internally corrects for library size, so transformed or
normalized values such as counts scaled by library size should not be
used as input.

Normalized counts like FPKM are not supported and should be converted to counts using the appropriate import funcion.
I don't really understand the use of the normalization factors for quality control because the number of reads in a sample does not correlate well with the quality. Maybe you could link to the paper in case there is something hidden in the experimental design that would explain this.
If you are interested in quality control a useful metric in my experience is the 3' sequencing bias that can be calculated using, for example, RNA-SeQC.


Error in if (is.na(sra_acc$run[i])) { : argument is of length zero

After checking with file.exists(), I would suggest you try the same command with another SRA run. It may happen that some runs are not accessible for reasons I do not know. I always have problems in accessing sra formatted files and use a workaround to download fastq.
Here is my not official workaround:
runlist<-sraConvert("SRP045534",sra_con=sra_con)
runlist[1:2,]

study submission    sample experiment        run
1 SRP045534  SRA179382 SRS685546  SRX681992 SRR1552451
2 SRP045534  SRA179382 SRS685544  SRX681990 SRR1552449

for (dloadrun in runlist$run)
{
    getSRAfile(dloadrun,sra_con=sra_con,fileType="fastq")
}

This should iteratively download all the runs associated to the project you are studying.
It is possible that you can do the same with sra type file, but I have never been able to use them with sraDB, and never investigated the issue.


Can blat use more than one core/CPU to speed up the alignment?

BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment.

EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based pipeline. So far there are no converters for minimap2's SAM/PAF or MUMmer's delta formats. With the UCSC pipeline, you have to write a converter by yourself. An alternative approach is to use minimap2's liftover. The minimap2 cookbook gives an example:

minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -

Note that this pipeline is not as mature as UCSC liftOver and is not as scalable.


What kind of analysis can be done with differential expression of transcription factors?

Expanding my comment.

TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. Perhaps there are some databases for TF.
I found a website base on the ENCODE project with the genes regulated by each TF.You can use them for pathway or gene set enrichment analysis, but I am not aware of any specific analysis for genes regulated by a single TF.

The histone data is harder if it was methylation I would say to look if different states of the histones are in different spatial state, but with gene expression is harder to conclude anything. Nevertheless you can check methylation studies to see if gene expression is associated with methylation changes.


Use Htslib to extract specific region of vcf file

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


Wow Biobash, that's kinda harsh. It is cool that sometimes you unexpectedly learned something new. What Michael mentioned is pretty interested, no one has ever mentioned that or I have never heard that.
So, my approach would be that you now are doing a Poolseq study, since you can't identify individuals. The other way to try to catch the paralogues is to map your reads to a sister species. My question is what are you trying to estimate? Genetic differentation? Genetic diversity? Pop Gen structure?
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6116966/
https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.5240
https://www.g3journal.org/content/9/12/4159
One semi-robust approach to assess usefulness is thinking about your system and defining what you expect to find, let's say in terms of genetic differentiation or genetic structure, and estimating that from your data. If that hypothesis is confirmed, your dataset might be useful. You can test the same by first estimating the loci by mapping your reads to a sister reference genome, to try to detect paralogues. Since you don't have barcodes, then each "locus" is independent, and you lose power, but there is nothing else you can do.
Good luck;


Using signal peptide and the expression levels of signal recognition particle in secretome analysis

This is something that Qiagen's Ingenuity Pathway Analysis is meant to solve. Unfortunately, the process of solving pathway dependency problems like these is difficult, and frequently leads to the discovery of super-linked proteins that are seen everywhere (and consequently fairly uninteresting from a biological perspective). The challenge is in filtering out the useless noise of "statistical significance", and identifying the biologically-significant signal. That's not something that computers are particularly good at doing at the moment.

My own suggestion is that you would be better off finding someone who has domain knowledge in the area of secreted proteins, and trying to present results to them in such a way that it makes it easy to interpret the results and ask further questions. If you are that person, perhaps there's a more specific question that you could ask to further investigate the problem that you have.


Extracting genes differentially expressed by Wilcox test

All of the genes are behind the difference, you're doing a rank test comparing all of them against each other at once. Your results have more to do with samples having slightly different outlier genes than any interesting biology. You don't have enough replicates to do a non-parametric test.


How to select specific lines and specific columns

You can simply use subset function in R and use the code given below,
    new_dataframe <-
    subset(name_of_your_dataframe,Column1 == 1 &&
    Column2 ==  "T1")

Please keep in mind that R is case sensitive and you should be careful while using the variable names.


Median string problem & multiple sequence alignment

Why is this approach not used? Are there any algorithms that use it?

As shown in this paper, the MEDIAN STRING problem is NP-Complete on alphabets of carnality greater than 3. While the problem isn't used directly in approaches due to the complexity, the idea of trying to create a consensus sequence between pairs of input and using that as building blocks is common. These are considered "progressive alignment methods".


Problems replicating Kover paper results

The single quote at 'mycobacterium is not matched so the command prompt is waiting indefinitely for the command to be completed. Maybe you want 'mycobacterium tuberculosis/kanamycin/metadata.tsv' (note closing quote)


I need some tips and suggestions for further analysis of NGS expression data (log2cpm)

Since you indicate already having performed differential expression and using gene set enrichment the next step is to look at what changes due to drug treatment in the pathways of interest while having a decent understanding of those pathways. Things like IPA or pathview (from Bioconductor) can aid in this, but at the end of the day you need to have a decent understanding of the pathways yourself to judge what seems to be like the best targets for future work.


If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.


Number of sing-copy genes identified by orthofinder is far less than BUSCO assessment?

It's hard to say for sure without seeing the exact commands you used. However, Orthofinder identifies/defines single-copy orthologues as those genes present in a single copy in all species being compared in the analysis. If a BUSCO orthologue is present in more than one copy in one of the species being compared (could be multi-copy, could be fragmented etc.), then it won't be recovered as single copy by Orthofinder.


Calculating the charge of a peptide computationally

A quick google search turns up protcalc, which is able to give a nice pH-dependent table of peptide charges (yours ranges from 3.1 at pH 4 to 1.5 at pH 10). It's on sourceforge, so hopefully the source code (and maybe a publication) is somewhere in there. Granted, this isn't R, but it's not clear from your post how necessary that really is.

BTW, in general "best and most accurate" turns out to be highly subjective (at least the "best" part of that).


Yes, this exists, and can be efficiently generated by plink 2.0's --export ind-major-bed command.  (The third byte is 0 instead of 1, and the specification is otherwise identical to that of regular plink .bed files with samples and variants swapped.)


Details of DESeq2 modeling a batch effect

DESeq2 uses the batch information (and everything else in the design) to produce offsets for its GLM. For a background on that please check how linear models work, e.g. using the StatQuest series of statistics videos over at YouTube.
It still operates on the raw counts. The same goes for the normalization factors.
removeBatchEffect fits a linear model to the data including the batch information and then subtracts the batch component from the counts (that is basically the baseline difference).
If you are interested in preserving the integer nature of counts and preserving zeros as smallest values after explicit batch correction you may want to check ComBat-seq() from the sva package. It operates on raw counts and returns batch-corrected raw counts which you could then normalize and calculate CPMs from (if you need batch corrected CPMs with no negative values). I find it useful and prefer it over removeBatchEffect as it avoids the unfortunate negative counts which sometimes messes up plotting scripts that expect zero as the smallest possible value.


Why doesn't the mouse GRCm38/mm10 refGene genome annotation file contain non-coding transcripts?

It wasn't the correct file.
I downloaded the file again from UCSC, compared it with the file that I was using, and they were different. This new file contained both information from coding and non-coding transcripts.
To all the kind people that took interest in my problem, I am sorry.
Edited 30/12/2020:
To add a bit of clarity; the link in the question is the correct link from which to download the correct file. The file with which I was working from the beginning was the wrong one, but I don't really know where I got it from the first time.


What do you want to do with NA values?

The below scripts ignore them (h3ab74's answer will print them to positiveExpression)

awk '$2 ~ /^[0-9]+$/  {print $0}' expressionFile.txt > positiveExpression.txt
awk '$2 ~ /^-[0-9]+$/ {print $0}' expressionFile.txt > negativeExpression.txt

These do a regex match to respectively, positive and negative numbers. If there's anything else in field 2 the line won't be printed.


the variation between treatments is less than the variation between replicates in RNA-seq data

If you got DEGs with statistically significant p-values after multiple testing following a valid pipeline, I'd believe that, even if looking at the first few PCs doesn't look promising.


Is it possible for coronavirus or SARS to be synthetic?

The scenarios are impossible and would be laughable if they were not so serious. The evidence is in the phylogenetic trees. Its a bit like a crime scene when the forensics team investigate. We've done enough crime-scenes often going to the site, collecting the pathogen, sequencing and then analysis - (usually neglected diseases) without any associated conspiracy theories.

The key technical issue is coronaviruses are zoonoses , pathogens spread to humans from animal reservoirs and phylogenetic tree really helps understand the  how the virus is transmitted.

Trees

The key thing about all the trees are bats. Bats lineages are present at every single point of the betacoronavirus phylogeny (tree), both as paraphyletic and monophyletic lineages, one example is this tree of betacoronaviruses here. Meaning the nodes connecting the branches of the tree to the "master-branch", represent common ancestors and these were almost certainly bat-borne cornaviruses. This is especially true for SARS and - here bat-viruses are EVERYWHERE.
The tree here also shows that SARS arose on independently on two occassions, again surrounded by bat lineages and 2019-nCov has emerged separately at least once, again associated with bats.
Finally the tree below is a figure from BioRxiv Zhou et al (2020) "Discovery of a novel coronavirus associated with the recent pneumonia outbreak in humans and its potential bat origin" shows the 2019-nCov lineage is a direct descendent of a very closely virus isolated from a bat (RaTG13*). This is a really conclusive finding BTW.

Note, I don't normally present inline images, but it is such a nice finding (hint to reviewers) and BioRxiv is open access.

Conspiracy theory 1: laboratory made virus

Literally it would require someone passaging a new virus, with unknown human pathogenicity, and independently introducing all the earlier passages enmass across bat populations of China. They would then hope each lineage becomes an indepedent virus population before, then introducing the virus to humans. Thus when field teams of scientists go around using mist nets to trap the bats, buy them from markets, isolate the virus and sequence it they would find a beautiful, array of natural variation in the bat populations leading up to the human epidemics, that perfectly matches vast numbers of other viral zoonoses. Moreover, this would have to have happen substancially prior to SARS and 2019-nCov, because the bat betacoronaviruses have been known about prior both epidemics, viz. its simply not feasible.

Biological explanation
General Bats are a reservoir host to vast numbers of pathogens, particularly viruses, including many alphaviruses, flaviviruses, rabies virus and beleived to be important in ebolavirus (I don't know about this) and even important to several eukaryotic parasites. It makes sense, they are mammals, so evolutionary much closer to us than birds for example, with large dispersal potential and roost in 'overcrowded' areas enable rapid transmission between bats.

Technical
The trees show bats are the common ancestor of betacoronaviruses in particular for the lineage leading into the emergence of 2019-nCov and SARS, this is seen in this tree, this one and the tree above. The obvious explanation is the virus circulates endemically in bats and has jumped into humans. For SARS the intermediate host, or possible "vector" was civet cats.

The theory and the observations fit into a seamless biological answer.

Conspiracy theory 2: Middle Eastern connection

I heard a very weird conspiracy theory attempting to connect MERS with 2019-nCov. The theory was elaborate and I don't think it is productive to describe here.

Biological explanation
All the trees of betacoronaviruses show MERS was one of the earliest viruses to diverge and is very distant from 2019-nCov, to the extent the theory is completely implausible. The homology between these viruses is 50%, so its either MERS or 2019-nCov. Its more extreme than mixing up yellow fever virus (mortality 40-80%) with West Nile virus (mortality <<0.1%), the two viruses are completely different at every level.

What about errors? Phylogeneticists can spot it a mile off. There are tell-tale phylogenetic signatures we pick up, but also we do this to assess 'rare' genetic phenomina. There is nothing 'rare' about the coronaviruses. The only anomaly is variation in the poly-A tail and that is the natural variation from in vitro time-series experiments. Basically we've looked at enough virses/parasites through trees, that have no conspiracy theories at all (often neglected diseases), and understand how natural variation operates - so a phylogenecist can shift the wheat from the chaff without really thinking about it.

Opinion
The conspiracy theories are deeply misplaced, and the only connection I can imagine is its China. However, the Chinese have loads of viruses, influenza in particular which causes major pandemics, but that is a consequence of their natural ecology (small-holder farming) allowing the virus to move between reservoir hosts. I've not visited small-hold farms in China, but I have in other parts of the world and when you see them, you get it. The pigs, chickens (ducks .. China), dogs, horses and humans all living within 10 meters of each other.

Conclusion
Shipping large numbers of bats to market, bat soup, raw meat from arboreal (tree-living) mammals such as civets that are sympatric to bats. Then consider the classical epidemiology in light of the phylogenetic data, which is very consistent, a single picture emerges that coronavirus is one of many zoonoses which has managed to transmit between patients.

Summary The fundamental point is the bioinformatics fit into the classical epidemiology of a zoonose.

*, Note The bat coronavirus RaTG13 predates the 2019-nCov outbreak by 7 years. It is not even clear whether the virus has been isolated, i.e. could just be a RNA sequence.

"They have found some 500 novel coronaviruses, about 50 of which fall relatively close to the SARS virus on the family tree, including RaTG13—it was fished out of a bat fecal sample they collected in 2013 from a cave in Moglang in Yunnan province."

Cohen, "Mining coronavirus genomes for clues to the outbreak’s origins" Feb, 2020 Science Magazine,


Appropriate tool or algorithm for sloppy alignment of degenerate bases

I believe mafft can align degenerate sequences.  From the documentation for version 7 at: https://mafft.cbrc.jp/alignment/software/anysymbol.html

Acceptable nucleotide symbols in the default mode:
a, c, g, t, u
Ambiguous nucleotides (r, y, w, s, k, m, d, v, h, b; IUPAC-IUB codes) can be used and are scored as:
score(r,a) = ( score(a,a) + score(g,a) ) / 2
score(r,g) = ( score(a,g) + score(g,g) ) / 2
score(r,t) = ( score(a,t) + score(g,t) ) / 2
score(r,c) = ( score(a,c) + score(g,c) ) / 2
...
score(r,y) = ( score(r,t) + score(r,c) ) / 2
...
score(r,r) = ( score(a,a) + score(g,g) ) / 2
score(y,y) = ( score(c,c) + score(t,t) ) / 2

MAFFT is designed for multiple sequence alignment, so will attempt to find a local alignment of homologous segments (using FFT), and then attempt to make a global alignment between the input sequences.


Yes. The general term for identifying subsets from a mixed group is deconvolution.

This is easier to do when there is haplotype-level information available (e.g. sequenced reads that span multiple SNPs), but it is possible to get some indication of sample mixes from unlinked SNPs alone.

Deconvolution is helped a lot by having proportional data in the input (e.g. frequencies of genotypes, rather than a single genotype). Most SNPchip platforms have multiple probes per variant, and can provide some indication of variant proportion when digging down into the raw data. Deconvolution is also a lot easier with some knowledge about the component mixes. It's much easier to simulate mixes than to guess at possible mixes.

Here's a paper about the application of deconvolution in forensics.

Here's a paper about a program called CIBERSORT which carries out a deconvolution of immune populations based on gene counts:

Here's another paper which discusses the idea in more general terms.

You might find more information on one of these Biostars posts:

Deconvolution Methods on RNA-Seq Data
Differential gene expession analysis in cell populations of mixed tumor and normal cells


This is a a symmetric matrix with the pathways in both axes. So for each pair a similarity index is calculated (usually it takes into account the number of genes in each pathway and the total number of genes in both pathways). There are several similarity indexes available. A typical one is Jaccard or Sørensen-Dice index, usually 0 means no similarity and 1 equal pathways.

You can do this kind of analysis in R using BioCor (Disclosure: I am the author of this package) which uses the Dice similarity. You will need the name of the genes of each pathway, after calculating the similarity of pathways using mpathSim, you can then plot it with a heatmap (using the heatmap function).


If I understand your question correctly, then I think for case of pairwise alignment, there is a simple explanation.

I believe the key insight is that: a mismatch should always score better than a gap.*

This follows biologically since the insertion/deletion (indel) rate is roughly 1/10th that of the substitution rate (i.e. the occurrence of single nucleotide changes), at least in vertebrates. (This varies across the tree of life but I think the substitution rate virtually always exceeds the indel rate.)

To understand why this matters, consider an example:

ATG-AG
ATGT-G

This is an 'impossible alignment' under the probabilities you gave since here we have a transition from a gap-residue alignment to a residue-gap.

However, under our assumption that mismatches are more likely biologically than indels, the correct alignment should be:

ATGAG
ATGTG

Indeed, the latter does look like a better alignment.

This also follows for more complex examples, so this:

ATG--AAG
ATGTT-AG

Becomes this:

ATG-AAG
ATGTTAG

(Or this:

ATGA-AG
ATGTTAG

)

* Strictly, I mean a substitution should score better than an indel (with the associated gap opening and extension penalties). In fact, for the assumption to always be true, a run of mismatches should still score worse than a single indel. This may not always be a correct assumption, consider this example below, is the true alignment case 1) or 2) or something else? Or is in fact a global alignment bad here and this should be split into 2 local alignments? Is there a likely biological mutational event that could explain this? I ask these questions just to point out it is not black-and-white, I don't have clear answers

1)

CGTACGTAGAGGAATGCCCCCCCCC--------AGCAACGTAGCAT
CGTACGTAGAGGAATG---------TTTTTTTTAGCAACGTAGCAT

2)

CGTACGTAGAGGAATGCCCCCCCCCAGCAACGTAGCAT
CGTACGTAGAGGAATGTTTTTTTT-AGCAACGTAGCAT


Collapse cell barcodes distribution within 1 Hamming distance

This doesn't use standard unix tools, but rather the API for UMI tools:

#!/usr/bin/env python
from umi_tools import network

cb = []  # cell barcodes
cbc = dict()  # cell barcode counts
for line in open("UMIs.txt"):
    cols = line.strip().split("\t")
    cols[1] = bytes(cols[1], "ascii")
    cb.append(cols[1])
    cbc[cols[1]] = int(cols[0])

uc = network.UMIClusterer()
CBclusters = uc(cb, cbc, 1)
cbFinal = dict()
for l in CBclusters:  # This is a list with the first element the dominant barcode
    cbFinal[l[0]] = 0
    for x in l:  # Iterate over all barcodes in a cluster
        cbFinal[l[0]] += cbc[x]
print(cbFinal)

If you give it input in UMIs.txt of the following form (note the last two lines differ from the first with a distance of 1):

11612552    TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
20  TCCAGAGCACTGCATAACTCAA
2   TCCTGAGCACTGCATAACTCNA

Then it will produce:

{b'TAAGGCGATCTCTCCGTACGTG': 5985219,
 b'TAAGGCGAACTGCATAACATGT': 6626630,
 b'GCTACGCTACTGCATAAGTCCA': 9349711,
 b'CGTACTAGTCTCTCCGCGGCTA': 8161950,
 b'TCCTGAGCGTAAGGAGCAGATC': 8102383,
 b'CAGAGAGGCTAAGCCTGCACAT': 8343678,
 b'TCCTGAGCACTGCATAACTCAA': 11612574,
 b'AGGCAGAATCTCTCCGAGTGCT': 6210446,
 b'AAGAGGCAGTAAGGAGATAATC': 7110298,
 b'CGAGGCTGTATCCTCTACTAGA': 6390489}


Comparing phylogenetic models with different datasets

The key issue is the data is subject to reticulation, languages are seldom clean breaks. In nucleotide data this can be identified due to deviations from physical linkage and call it recombination. Linkage doesn't apply to languages.

The fundamental assumption is one tree - the same tree - describes all the data, with reticulation that may not be a good assumption. Thus English is a hybrid language, and lots of English words find their way back into French and German 'parents'. Or perhaps a better example is Gaelic mingles with English over time. There are methods used to monitor reticulation in languages.

Thus id keep it simple. The branch length is not only directly comparable it provides linguistical information that is useful to understand. This is the amount of 'evolutionaey' change per word (it that is what you binary data is describing) and the relationship to the node ('breakpoiint') is of primary interest.

The second step is to simply draw a tree of data set A and compare it to data set A+B and data set B and the probabilities describing the sister groups. Are the trees congruent, i.e. are they statistically the same?

The area you describing is more given to parametric bootstrapping, using custom resampling, but you need to be certain that the model being used is appropriate for the data.

The final point is whether a reversible or directional matrix is being used, which impacts on the specific result.


16S species-level taxonomic assignment--what is the current state of the art?

I have been using dada2 and phyloseq (phangorn for trees) combination to analyze a V4 region 16s rRNA data set. I find their approach and explanations very understandable and usable since I prefer R.
I was under the impression that even full-length 16s rRNA gene sequences were inadequate to classify up to species level since closely related species within a genus can have identical 16S rRNA gene sequences.
But, turns out, intra-genomic variants of this gene within a species or even strain may allow their unambiguous identification according to this paper:
Johnson, J.S., Spakowicz, D.J., Hong, B. et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun 10, 5029 (2019). https://doi.org/10.1038/s41467-019-13036-1) available at: https://rdcu.be/b5cgB
Another relevant article regarding using black box ML algorithms (random forests) versus a straightforward logistic regression is discussed in this one:
Topçuoğlu BD, Lesniak NA, Ruffin MT 4th, Wiens J, Schloss PD. A Framework for Effective Application of Machine Learning to Microbiome-Based Classification Problems. mBio. 2020;11(3):e00434-20. Published 2020 Jun 9. doi:10.1128/mBio.00434-20
Since the OP has probably already moved on long since the original query, I am writing this to document relevant current references for myself and anyone researching this topic on StackExchange Bioinformatics.


What is the purpose of 'Folder locations' in MaxQuant?

These folders are usually generated in a default directory (where the raw files are). It can however be beneficial to set the paths to these folders manually.

  Temporary folder

that is obviously where temporary files go, so having an SSD or other fast drive will benefit the I/O operations. Set that path to a fast storage partition.

  Combined folder

all the output, including txt directory will go here. Setting it to fast storage might speed up writing the txt files. If you set it to the same path while running multiple instances of MQ it will obviously override output.

  Andromeda index

this is where index files from Andromeda search are. By default these are inside a combined/search folder. Since these index files can get big, once generated by MQ they can be reused. Set this path if you want to reuse previously generated index files.


Get RefSeq accession numbers with versions

I don't believe this is possible using biomaRt, nor using AnnotationHub.

I have two suggestions, neither of them very satisfactory. First, you can specify an Ensembl archive for biomaRt, for example:

mart72.hs <- useMart("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl",
                      host = "jun2013.archive.ensembl.org")

Of course, that requires that you have some idea of the date for each accession version and that the archives span that date - so not especially useful.

The other option is to access EUtils using e.g. rentrez, which does allow search by version number:

library(rentrez)
es <- entrez_search("nuccore", "NM_002084.1 NM_002084.2")
es$ids

[1] "4504104" "6006000"

So knowing the accession, you could simply append 1, 2, 3... to it, run the search and see if UIDs come back, then get them using entrez_fetch.

EDIT: here's a quick and dirty function which takes an accession as query, appends version = 1, fetches the ID, then increments the version and repeats until no more results are returned. It is not well-tested!

getVersions <- function(accession) {
  require(rentrez)
  ids <- character()
  version <- 1
  repeat({
    es <- entrez_search("nuccore", paste0(accession, ".", version))
    if(length(es$ids) == 0) {
      break
    }
    ids[version] <- es$ids
    version <- version + 1
  })
  ids
}

Example:

getVersions("NM_002084")
[1] "4504104"    "6006000"    "89903006"   "1048339180"


Major variability in speed of BLAST between proteins

Try changing you evalue threshold to something like 0.001.

If you blasted a sequence against a database containing 10,000 sequences and got a hit with an evalue of 1, this would mean that, given the size of your database and length of protein, you would expect to find a match as good as (or better than) the hit returned once purely by chance. If the evalue threshold in your query was set to 0.01 (for example) then this hit would not be returned.

If you got a hit with an evalue of 200,000 this would mean you would expect to see 200,000 matches with a similar score simply by chance. If your evalue threshold were set to 200000, then all of the hits with an evalue <=200000 would be considered significant. If you set it to 0.001 then only those hits with evalues <= 0.001 would be considered significant.

This means that (depending on the size of the data base) almost every sequence in your query set will be considered as significant given your threshold, which will slow down your search hugely.

See the the section on E values here


Documentation for Transfrag class codes (cuffcompare)

You can get a bit more information by looking at the source code comments. Contained is for single exon transcripts that are primarily contained (>80% overlap) within an exon of an already annotated transcript. Generic exonic overlap is for cases where a transfrag exon overlaps that of an annotated transcript. The strand doesn't seem to matter in such cases.


What score value is considered acceptable when using Vina Autodock tool for molecular docking simulations?

The answer is "Depends". But the best thing is always to look at what you ended up with and do some triaging before buying/synthesising anything. By triaging I mean filtering by its properties (.e.g. Lipinski, reactivity, synthetic feasibility etc. —discussed in this Bioinfo SE question).
Score
Regardless of docking software used, a forcefield free-energy score is a result of the scoring function used and of the properties of the molecule —including VENA which uses a hybrid score function.
This score is not, but generally is set up to try to match, the actual phys-chem one.
Gibbs free energy is measured in kcal/mol (or less commonly in kJ/mol). A hydrogen bond is about 2 kcal/mol, a π-π bond 1 kcal/mol and a salt bridge 3. The average energy of a collision of RT water molecules is 0.6 kcal/mol.
Therefore, if you dock multiple compounds (virtual library screen, VLS) you sort by free energy to see the strongest binders and consider the top ones.
However, a catch is that this will be biased by MW.
Also, it is useful to score the hits with a different program (say Gold or ICM-dock).
Ligand Efficiency
A variant metric that is employed when your hits vary in size and when they are over a given size (e.g. not fragments under 200 Da) is Ligand Efficiency (LE), which is the ratio of your score over the number of heavy atoms (ref).
Single compound
If you have a single compound and want to know if it binds acceptably, then this is a length-of-a-piece-of-string question as every compound has an affinity constant to a given protein, except most may be beyond their solubility. Initial hits will have µM affinity, while a lead will have nM affinity. However, a KD is not a ∆G and depends on your protein.


How to identify the reference genome assembly used to create a .hic file?

Turns out HiC file does store such metadata.

straw tool from Aiden Lab contains a script (read_hic_header.py) that prints out the metadata from a given HiC file. Usage: read_hic_header.py <hic file or URL> [verbose]

In order to obtain the metadata in the form of a dictionary, one may consider using a function of the same tool that is currently available in my fork.

from straw import straw_module
straw_module.read_metadata(infile,verbose=False)

Args
    infile: str, path to the HiC file
    verbose: bool, print the metadata if True else return the metadata
Returns
    metadata: dict, the metadata is returned if verbose is False.
                Keys of the metadata:
                HiC version,
                Master index,
                Genome ID (str),
                Attribute dictionary (dict),
                Chromosomes (dict),
                chromosomes2length (dict),
                Base pair-delimited resolutions (list),
                Fragment-delimited resolutions (list).


What are the optimal parameters for docking a large ligand using Hex?

While I have no experience with this specific question you have, according to Feinstein & Brylinski (2015), a "fully automated procedure", i.e. a Perl script, can be used to optimize the box size itself, and it can be found here.


I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.

No, a phenotype is a way of behavior or other observed characteristic of a person, resulting from a combination on its underlying biology and the environment it's in. Think things like, height, running speed, strength, etc.

In the context in which you're interested, this would be more properly, "the clinical presentation of a disease or disorder".

Step one of diagnosing a patient is defining exactly how they present (i.e., exactly what their maladies, aberrant electrolyte levels, genetic mutations, etc. are).
An algorithm that only works on the test population is useless. Google "model overfitting".
The point of an algorithm is to take the EHR data and progress to a diagnosis. This isn't real time, you'd just run a batch of algorithms as test results come in.
You need patients in your test/validation sets to have the required inputs for your algorithm.

FYI, I don't know of any regulars on this site that deal with EHR data. If there's a community of such people on the web they must be elsewhere.


Why are the clusters in the left side of a tSNE plot and not in the centre?

I'm going to stick my neck out here. Overall PCA followed by tSNE is a very powerful method of unsupervised learning however, I agree with @GWW,

cluster <- ScaleData(object = cluster, display.progress = FALSE)

If this is a max-min rescale transformation (common in supervised learning) then the outliers will badly affect the clustering. You can remove outliers that are over 2 sd from a formal analysis and re-run the calculation. That is just the risk of a max-min rescaling.

The one issue I'm concerned about is the initial quasi-standardisation followed by log-transformation,

log1p(value/colSums[cell-idx] *scale_factor)

In this case you are multiplying by 10000, this is to remove the ratio (used as a quasi-standardisation) and then taking the log ??? That is not the way to transform a ratio, it would be better to take the reciprocal of the ratio. There are other (better) methods, a formal standardisation would be useful, before a transformation is considered. Possibly if that was done correctly you wouldn't have the problem of outliers potentially crunching your data.

If you fiddle with the parameters you might get something - but its a symptom not the cause of the problem. However, if removing the outliers works (I have not checked the ScaleData function), then you are sorted and all else is forgotten (at least for this data set). The information I read was here.


Importance of Proper Pairs vs Aligned Reads for RNASeq data

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%).


How to generate two fasta files from samtools phase output?

Answer January 20th, 2020. If you run samtools phase with the option -b my_prefix it outputs a bam file for both haplotypes, named my_prefix.0.bam and my_prefix.1.bam. If the sequence has been correctly phased, use these bam files to correct the reference using pilon. This is preferable to the method described in the January 14th, 2020 answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that samtools phase found. samtools phase does not find indels.

Answer January 14th, 2020. This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.

Run the program like this:

samtools phase reads.bam | phase2fasta.py ref.fa > haplotypes.fasta.

Here is the code of phase2fasta.py:

#!/usr/bin/env python

import os
import sys

if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
    raise Exception("The first and only argument is the reference fasta.")
allele0 = []
seqname = ""
with open(sys.argv[1]) as f:
    for line in f:
        if line[0] == ">":
            seqname = line[1:].split(" ")[0].strip()
        else:
            for c in line.strip().upper():
                allele0.append(c)
allele1 = [x for x in allele0]
for line in sys.stdin:
    if line[0:2] == "M1":
        fields = line.strip().split()
        if fields[1] == seqname:
            pos = int(fields[3]) - 1
            allele0[pos] = fields[4]
            allele1[pos] = fields[5]
print(">{}_allele0".format(seqname))
print("".join(allele0))
print(">{}_allele1".format(seqname))
print("".join(allele1))


How to visualize called narrowPeak files in UCSC Genome browser?

Headers should not contain tabs:
track type=narrowPeak name="Somite narrowPeak" description="Somite narrowPeak"

Ensure that you have NO tabs on that line.


Metagenomic shotgun data with internal control

The problem with correlations in relative scale is know on the literature, it has been suggested to change to other metrics. So I wouldn't use it for internal control.

As you can see on figure 3 of the linked article the correlation on the absolute scale tend to be higher. So I would say this is within expectations.

As other way for internal controls are blanks (Just water to see if there are any contaminants), and technical replications, to make sure that the measure is correct not some problem with picking some amplifications instead of others...

Also you don't mention it but the tool and library used to impute the microorganism is also important. You can check with other tools to see if the estimation is similar.


Remove repetitive region vcf file using repeatmasker bed file

A couple fast alternatives:
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(sort-bed repeats.bed) > answer.bed
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(rmsk2bed < repeats.rmsk) > answer.bed

Using --chrom chr14 focuses set operations on that chromosome.


should I continue the iteration of psiblast if the blast step produced only below threshold sequences?

The first round of psiblast is the same as blastp and an evalue of 0.001 is low and below the recommended threshold of (e)-value = 0.005.

You require significant knowledge of composition bias and globular domains to move into the 2nd and 3rd iterations. I'd recommend reading NCBI's manual here.

From you response you are using psiblast as a tool to look at antibody cross-reactivity, which I presume is pan-Clostridium.

I  recommend using psi-blast as investigative analysis only and switching to the Immune Epitope Database here. Everything will be explained on the IEDB server. On this server I suggest performing two separate analyses and seeing whether the results correlate,

B-cell epitope prediction on your strongest evalue
B-cell epitope prediction on your query sequence
Do the epitopes positions match?

If the epitope positions do not match then the hypothesis of a pan-Clostridium cross-reactive antibody can be immediately discarded. If there is a match it provides the basis for further investigation.

Please keep in mind Clostridium will use glycosylation to modulate antibody binding.


server for finding kmers in set of sequences

One Codex is powerful and user friendly, and works in a similar way to Kraken underneath
https://www.onecodex.com/

The Kaiju web server might be of interest and is more sensitive than exclusively k-mer based methods
http://kaiju.binf.ku.dk/server

The title of this question is misleading – it sounds like you're actually looking for a webserver for taxonomic classification using k-mers.


How can I dock a protein to a nucleic acid?

I'm not sure what you meant but you can take a look at NPDock (disclaimer we wrote that tool). If you have a structure of your protein of interest, you can dock it to the structure of your DNA/RNA of interest. Mind that this is a rigid body docking which means that the structure will not change upon binding.


Is there a modern alignment tool tailored for transmembrane regions?

Answer from @matteo-ferla converted from comment:
I assume you mean sequence (MSA), not structure alignment. Despite being a structural bioinformatician, I would slap a big caveat against structure based aligners. TCoffee:Espresso is probably the best structure-based aligner (browser and download and API), but the results are not as solid as some non-structural approaches such as training Psi-blast and aligning with that matrix. I am not familiar with MSA of TM regions, however, so this group may be an exception.
Given that helices are a repeat of 2 hydrophobic residues and 1.5 either polar (core-wards) or also hydrophobic (depending on how the helices pack), a pattern based aligner ought to do well from an evolutionary perspective. If you really want from a structure perspective and don't want to know about the cases where a loss of a turn at the top is compensated by the bottom, say, then TCoffee is definitely worth a try
Answer from @m converted from comment:
PSI-blast is a useful tool for deep alignment patterns, albeit it has its critics. I agree the vast majority of in vivo structures lack TM proteins, for different reasons


I don't know if this is the best answer, but one quick-n-dirty approach I commonly use involves tweaking bwa mem to report only perfect matches. So in your case, the idea is to make two copies of each 21-mer, replace the N with the appropriate allele, and then map the resulting sequences to your bwa-indexed assembly like so.

bwa mem -k 21 -T 21 -a -c 5000 assembly.fasta your-21mers.fasta

The -k 21 says to use a minimum seed length of 21, which forces exact matches. The -T 21 requires a minimum score of 21, which also enforces an exact match. The -a parameter reports all matches, since a "best" match doesn't make sense in this situation. The -c parameter limits how many matches are reported, which may need to be adjusted depending on how repetitive the 21-mer is.

UPDATE If you don't want to introduce a dependency like bwapy, it's fairly straightforward to implement this in pure Python.

with tempfile.TemporaryFile() as kmerfile:
    for kmer in kmers:
        print('>{k}\n{k}'.format(k=kmer), file=kmerfile)
    kmerfile.flush()             # Often works without these two lines
    os.fsync(kmerfile.fileno())  # But best to be safe
    bwacmd = ['bwa', 'mem', '-k', ksize, '-t', ksize, '-a', '-c', '5000', refrfile, kmerfile.name]
    bwaproc = subprocess.Popen(bwacmd, stdout=subprocess.PIPE, universal_newlines=True)
    for line in bwaproc.stdout:
        # process SAM records line by line

Also, one can call bwa mem via a C API, but the mechanics of this are left as an exercise for the reader.


Database with antigens and current vaccines

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.


How to model the effect of acid/base on pH

I figured that one has to include only the changes in the rules formula. Keeping the initial value of the Species or Parameter when defining circular or self-referencing formulas caused the error. So the correct rule becomes:

pH=-a*Acid instead of pH=pH-a*Acid or pH=7.4-a*Acid.

For convenience, I include the correct code crowded with several reactions producing both acids and bases, where pH gradually decreases as expected. The names are fictional, and just are better placeholders instead of x1 and x2, etc.

library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Pyruvate", 0)
addSpecies(glu, "Oxalate", 0)
addSpecies(glu, "Acid", 0)
addSpecies(glu, "Base", 0)
addSpecies(glu, "pH", 7.4)

addMAreaction(glu, "Glucose = Lactate + Acid", "forwardL",
              "backwardL", name = "glycolysis")
addMAreaction(glu, "Lactate = Pyruvate + Acid", "forwardP",
              "backwardP", name = "lactolysis")
addMAreaction(glu, "Pyruvate = Oxalate + Base", "forwardO",
              "backwardO", name = "oxalate synthesis")

addMAreactRate(glu, "forwardL", "assigned",
               "k_glucose_to_lactate*Glucose")
addMAreactRate(glu, "backwardL", "assigned",
               "k_acid*pH+k_lactate_to_glucose*Lactate")
addMAreactRate(glu, "forwardP", "assigned",
               "k_lactate_to_pyruvate*Lactate")
addMAreactRate(glu, "backwardP", "assigned",
               "k_acid*pH+k_pyruvate_to_lactate*Pyruvate")
addMAreactRate(glu, "forwardO", "assigned",
               "k_pyruvate_to_oxalate*Pyruvate")
addMAreactRate(glu, "backwardO", "assigned",
               "k_base*pH+k_oxalate_to_pyruvate*Oxalate")

addParameters(glu, "k_glucose_to_lactate", 0.02)
addParameters(glu, "k_lactate_to_glucose", 0.001)
addParameters(glu, "k_lactate_to_pyruvate", 0.01)
addParameters(glu, "k_pyruvate_to_lactate", 0.0001)
addParameters(glu, "k_pyruvate_to_oxalate", 0.01)
addParameters(glu, "k_oxalate_to_pyruvate", 0.001)
addParameters(glu, "k_acid", 0.0001)
addParameters(glu, "k_base", 0.0001)

addRule(glu, "pH rule", "ODEs", "pH=-k_acid*Acid+k_base*Base")

makeModel(glu)
simResults <- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")

The head of the resulting data.frame is:

  time   Glucose  Lactate   Pyruvate      Oxalate     Acid         Base       pH
1 0.00 1000.0000   0.0000  0.0000000 0.000000e+00   0.0000 0.000000e+00 7.400000
2 0.01  844.0001 155.0615  0.9384259 1.402282e-05 156.9384 1.402282e-05 7.399914
3 0.02  785.0869 210.3646  4.5476263 7.892761e-04 219.4615 7.892761e-04 7.399719
4 0.03  772.1557 218.6462  9.1924839 5.655403e-03 237.0425 5.655403e-03 7.399489
5 0.04  768.8212 217.2720 13.8876293 1.917807e-02 245.0856 1.917807e-02 7.399248
6 0.05  767.0383 214.4883 18.4278053 4.550630e-02 251.4350 4.550630e-02 7.399000

and resulting graph looks like:


Why do I obtain very strange results with DESeq2?

I'm guessing here, but I suspect this is at least in part due to the fact that all your samples in condition 1 are not expressed at all.

Formally you have an infinitely large fold-change between condition1 and condition2. Most tools solve this by adding a small pseudocount to unexpressed samples, but this will still result in a very high fold change. So far for the big fold-change.

If I had to guess at the reason why this is significant I would think in terms of dispersion estimation. This Bioconductor post has some discussion with the authors of DESeq2 on how this works.

Finally: calculating differential expression on these genes in these samples will never be really informative, as they really are only expressed in one sample. The way I see it you have two options:

Exclude the sample from your analysis. Straightforward, but does not solve this problem for other samples.
Filter your genes before calculating differential expression, i.e. remove all genes that are not expressed in a large number of samples. This approach is standard in for instance sleuth. A good starting point would be to remove genes that are not expressed in 75% of your samples.


Microbial Classification with 16S Nanopore Data

Use kraken2 with one of the 16S databases. Kraken2 has worked well for me with nanopore sequences. Pavian can be used to inspect the Kraken2 report files visually.

Using uncorrected called reads has worked for me. Any errors from nanopore sequencing should be inconsistent enough (i.e. not always lead towards the same wrong organism) that the noise doesn't matter.

Note: if you're using the MiniKraken database instead of a comprehensive 16S database, there might be some false positive species assignments depending on which kmers get pruned out.

Another possibility is centrifuge, but I'm not sure if that works for 16S amplicon sequencing; it's primarily for metagenomic analysis.


It depends of course on your available RAM on your computer.

Under Windows you can increase the available RAM in R with memory limit:

> memory.limit()
[1] 8070
> memory.limit(size=10000)
[1] 10000
> memory.limit()
[1] 10000

But that is not necessary on linux/unix.


Uppercase vs lowercase letters in reference genome

What does this soft masking actually mean?

A lot of the sequence in genomes are repetitive. Human genome, for example, has (at least) two-third repetitive elements.[1].

These repetitive elements are soft-masked by converting the upper case letters to lower case. An important use-case of these soft-masked bases will be in homology searches: An atatatatatat will tend to appear both in human and mouse genomes but is likely non-homologous.

  How confident can I be about the sequence in these regions?

As you can be about in non soft-masked based positions. Soft-masking is done after determining portions in the genome that are likely repetitive. There is no uncertainty whether a particular base is 'A' or 'G', just that it is part of a repeat and hence should be represented as an 'a'.

  What does a lowercase n represent?

UCSC uses Tandom Repeat Finder and RepeatMasker for soft-masking potential repeats. NCBI most likely uses TANTAN. 'N's represents no sequence information is available for that base. It being replaced by 'n' is likely an artifact of the repeat-masking software where it soft-masks an 'N' by an 'n' to indicate that portion of the genome is likely a repeat too.

[1] http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002384


I believe that

uniqsampleID = 1:length(unique(mutations$sampleID))
names(uniqsampleID) = unique(mutations$sampleID)
mutations$sampleID = uniqsampleID[mutations$sampleID]

should work, assuming mutations is a data frame.

Explanation:

1:N creates a vector of length N, containing 1, 2, ..., N. So using 1:length(X) will do this for each item in X, so 1:length(unique(mutations$sampleID)) should give you a vector of 1:139.

Then names(uniqsampleID) = unique(mutations$sampleID) assigns to each of these numbers one of these sample IDs, which allows you to retrieve the sample number with uniqsampleID[mutations$sampleID].


from https://www.blopig.com/blog/2019/08/how-to-iterate-in-pymol/
Example 1: Print the names of all the selected objects
You’ve selected several objects in PyMOL, and you want to see their names. By default, a new selection object called “(sele)” is created. We can create a Python variable called seleobjs, say, and then loop over them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: print obj

Example 2: Disable all the selected objects
If you still have your selection, “(sele)”, you can disable them (effectively turn them off on the right-hand sidebar):
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.disable(obj)

Example 3: Delete all the selected objects
You’re done with the objects in your selection, “(sele)”: you might be able to delete them using “delete ”, but if you can’t write a wildcard that covers all the names of the objects, you can iterate over a selection and delete them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.delete(obj)

from https://pymolwiki.org/index.php/Get_Area
get_area calculates the surface area in square Angstroms of the selection given. Note that the accessibility is assessed in the context of the object(s) that the selection is part of. So, to get the surface areas of e.g. a component of a complex, you should make a new object containing a copy of just that component and calculate its area.
please let us know if it worked


How can I use UpSet plot for expression/RNA-seq data

You need to generate a list of gene ID vectors to feed into upset. In the example below, I am selecting observations (genes) which have an level above the median level. You could also choose only genes above the 3rd quartile, etc. You then need to convert these from a list of character vectors to a presence/absence matrix. UpSetR has a fromList function that does just this.

In future, be more specific about what you want to achieve and supply example data for people to work with.

library("UpSetR")

tmp <- mtcars

tmp$car <- rownames(mtcars)
list <- lapply(colnames(mtcars), function(x) {
  tmp[tmp[[x]] > median(tmp[[x]]), "car"]
})
names(list) <- colnames(mtcars)

upset(fromList(list))

Created on 2019-10-17 by the reprex package (v0.3.0)


It seems to be some kind of a memory leak in pandas's read_csv, maybe due to it doing gzip-decompression.

Doing the decompression in a separate process and piping the result into cooler cload pairs should do the trick.

gzip -dc GSE109344_monkey_fibro_allValidPairs.txt.gz | cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} - /data/${FILE}.cool


How to extract the protein fasta file from a genbank file?

One can get it to work by using SeqIO.InsdcIO.GenBankCdsFeatureIterator:

from Bio import SeqIO

file_name = 'NC_000913.3.gb'

# stores all the CDS entries
all_entries = []

with open(file_name, 'r') as GBFile:

    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

This will then give the desired output (just the first two entries):

>b0001
MKRISTTITTTITITTGNGAG
>b0002
MRVLKFGGTSVANAERFLRVADILESNARQGQVATVLSAPAKITNHLVAMIEKTISGQDA
LPNISDAERIFAELLTGLAAAQPGFPLAQLKTFVDQEFAQIKHVLHGISLLGQCPDSINA
ALICRGEKMSIAIMAGVLEARGHNVTVIDPVEKLLAVGHYLESTVDIAESTRRIAASRIP
ADHMVLMAGFTAGNEKGELVVLGRNGSDYSAAVLAACLRADCCEIWTDVDGVYTCDPRQV
PDARLLKSMSYQEAMELSYFGAKVLHPRTITPIAQFQIPCLIKNTGNPQAPGTLIGASRD
EDELPVKGISNLNNMAMFSVSGPGMKGMVGMAARVFAAMSRARISVVLITQSSSEYSISF
CVPQSDCVRAERAMQEEFYLELKEGLLEPLAVTERLAIISVVGDGMRTLRGISAKFFAAL
ARANINIVAIAQGSSERSISVVVNNDDATTGVRVTHQMLFNTDQVIEVFVIGVGGVGGAL
LEQLKRQQSWLKNKHIDLRVCGVANSKALLTNVHGLNLENWQEELAQAKEPFNLGRLIRL
VKEYHLLNPVIVDCTSSQAVADQYADFLREGFHVVTPNKKANTSSMDYYHQLRYAAEKSR
RKFLYDTNVGAGLPVIENLQNLLNAGDELMKFSGILSGSLSYIFGKLDEGMSFSEATTLA
REMGYTEPDPRDDLSGMDVARKLLILARETGRELELADIEIEPVLPAEFNAEGDVAAFMA
NLSQLDDLFAARVAKARDEGKVLRYVGNIDEDGVCRVKIAEVDGNDPLFKVKNGENALAF
YSHYYQPLPLVLRGYGAGNDVTAAGVFADLLRTLSWKLGV


No counts for added gene in cellranger (scRNA-seq)

The first column of your GTF must correspond exactly to the name of the sequence in your FASTA file. Typically this is the chromosome so a GTF line:

chr1    hg38_knownCanonical exon    11106535    11262507    0.000000    .   .   gene_id "gene1"; transcript_id "tx1";

This must have a corresponding sequence in the FASTA file:

>chr1
ATCG ATCG ATCG ATCG ATCG ATCG...

Therefore you need to add an additional sequence to your FASTA file with >TdTomato. This must have exactly the same sequence as the construct that you are expecting reads for. The run cellranger mkref and cellranger count on these FASTA and GTF files. Cellranger also can take a GTF with “CDS” as a feature in addition to “exons”. The STAR alignment algorithm used by cellranger is splice aware so it can handle split reads but you will need to ensure the GTF corresponds to exons of the mRNA/cDNA rather than the full sequence.

10x Genomics is a platform with a high dropout rate, especially for lowly expressed genes. You should bear this in mind for what you expect to observe. It may be worth comparing the mapping rate of your construct to housekeeper genes to see if this is the reason for the low counts. Since whole mRNA molecules are captured on the beads before preparing the library, there should not be a difference in UMI counts between long and short genes.


Can you give a little more information about what tools you are using? How did you generate the OTU table etc?
FastTree is a very simple tool to generate a decent phylogeny. If you have the sequences of each OTU, you can align them in MegaX with muscle, or simply using muscle from the command line if you prefer. You can then generate the tree with this command (presuming you are on windows):
FastTree.exe -gtr -nt < aligned_sequences.fasta > otus_tree.tre

I hope that helps!


How to simulate NGS reads, controlling sequence coverage?

I am working on a Illumina sequencing simulator for metagenomics: InSilicoSeq

It is still in alpha release and very experimental, but given a multi-fasta and an abundance file, it will generate reads from your input genomes with different coverages.

From the documentation:

iss generate --genomes genomes.fasta --abundance abundance_file.txt \
    --model_file HiSeq2500 --output HiSeq_reads

Where:

# multi-fasta file
>genome_A
ATGC...
>genome_B
CCGT...
...

# abundance file (total abundance must be 1!)
genome_A    0.2
genome_B    0.4
...

I didn't design it to work with coverage but rather abundance of the genome in a metagenome, so you might have to do a tiny bit of math ;)


how to export the monocle data into Seurat

There's no default for the export_to option, so you just need to specify export_to="Seurat". Note that monocle2 doesn't seem to be compatible with the newest Seurat (3.1) releases.


using python to write bioinformatics pipelines tutorial

Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and  runs it robustly.

Snakemake is my favourite. It's (nearly) pure Python and can generate reports.
Nextflow is growing in popularity and is pretty straightforward
Ruffus used to be reasonably popular, seemed fine when I used it
Bpipe is for bioinformatics
Airflow is a more industrial "big" solution but in wide use
See this big list of pipeline systems

I'm sure you can find something there.


Whereas wouldn't be as fast as awk, here is a rather dirty R solution:

library(dplyr)
library(tidyr)
library(magrittr)

# File 1

# the "fill" parameter is used since not all the rows have the same number of cols
file_1 <- read.table("file_1.txt", skip = 1, fill = TRUE)

# paste() function with argument collapse is applied to the "phenotype"
# columns of each row (these cols start at position 6 and goes until
# the number of cols: dim(x)[2] gives the number of cols of a table
phenotypes <- apply(file_1[, 6:dim(file_1)[2]], 1, function(x) paste(x, collapse = ""))

# Cols referring to the genomic position, ref and alteration (cols 1:5)
# are combined with the newly created "phenotypes" column
file_1 <- data.frame(file_1[, 1:5], phenotypes)

# colnames are added, could be avoided if dealt with at file reading step
names(file_1) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# File 2

file_2 <- read.table("file_2.txt", skip = 1, fill = TRUE)

phenotypes <-apply(file_2[, 6:dim(file_2)[2]], 1, function(x) paste(x, collapse = ""))

file_2 <- data.frame(file_2[, 1:5], phenotypes)

names(file_2) <- c("Chr", "Start", "End", "Ref", "Alt", "Phenotypes")

# The two tables generated as above are merged using the appropriate
# columns as "anchors" in the "by" argument of full_join()
# The resulting Phenotypes.x and Phenotypes.y cols are concatenated
merged_tables <- full_join(file_1,
                           file_2,
                           by = c("Chr", "Start", "End", "Ref", "Alt"))
merged_tables$Phenotypes <- paste(merged_tables$Phenotypes.x,
                                  merged_tables$Phenotypes.y,
                                  sep = ",")
merged_tables$Phenotypes.x <- NULL
merged_tables$Phenotypes.y <- NULL

# Duplicated rows (in terms of genomic pos, ref and alteration) in any
# of the files result in duplicated rows in the merged table and are
# "handled" by aggregating/summarizing the info at the phenotype cols
# by using paste()
# strsplit() is used to "split" the "phenotypes" on "," in order to
# remove repeating phenotypes resulting from the full_join() call

merged_tables %<>% group_by(Chr, Start, End, Ref, Alt) %>%
      summarize(Phenotypes = paste(Phenotypes, collapse = ",")) %>%
      ungroup() %>%
      mutate(Phenotypes = strsplit(Phenotypes, ","))

merged_tables$Phenotypes <- lapply(merged_tables$Phenotypes, function(x) paste(unique(unlist(x)), collapse = ","))
merged_tables$Phenotypes <- unlist(merged_tables$Phenotypes)

# A tibble: 3 x 6
    Chr Start   End Ref   Alt   Phenotypes
  <int> <int> <int> <fct> <fct> <chr>
1     1  1000  1000 A     T     PhenoA,PhenoB
2     2  2000  2000 T     G     PhenoF,PhenoG,PhenoH,PhenoI
3     3  3000  3000 C     G     PhenoL,PhenoM,PhenoN

# The resulting table (tibble in this case) can be saved as a file
# using write.table() or write.xlsx()

paste() in combination with its collapse parameter is used because of the unequal numbers of columns in the files. fread() of the data.table package could be helpful for reading large files.


How to refine this distribution comparison script

I am not sure how useful this answer is. But may be you could try this.
For your python script:
You use f_out.write within a loop. I think usually this I/O operation is resource intensive. I guess the speed would improve significantly if you store the string in a variable and write the output only once outside the for loop like this:
out_str = ''
for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            out_str +=  (str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

f_out.write(out_str)

Regarding the memory error in your R script, are you using 32bit build of R?
Please see this:
http://web.mit.edu/~r/current/arch/i386_linux26/lib/R/library/base/html/Memory-limits.html#
They say "Under most 64-bit versions of Windows the limit for a 32-bit build of R is 4Gb: for the oldest ones it is 2Gb. The limit for a 64-bit build of R (imposed by the OS) is 8Tb."
Conclusion: Please try using your python script, but use the write operation outside the loop.


SNP-phenotype association analysis only using the SNPs on a specific gene

I know that the GWAS association p-value threshold is 1e-8

This may be a common threshold of statistical significance that is used, but it's definitely not an absolute value. It's a hack to try to work around many issues with GWAS associated with testing millions of SNPs. Unfortunately, the most relevant issue in GWAS (for spurious significance) is unexpected shared genetic structure in the cases or controls, and this cannot be excluded by a p-value threshold, regardless of how low it is set.

In your 100-SNP example, you are correct that the bonferroni correction would be a threshold of 1e-4 (for a desired significance threshold of 0.01), but it's important to be cautious about GWAS results, even when the p-value indicates otherwise. Trust no one, especially yourself. P-values should never be used as an indication of importance (which suggests that most GWAS results should probably be reconsidered), and are better used in conjunction with other evidence for identifying relevant / significant tests.

My recommendation for GWAS (or any similar bulk genotyping study, as in your 100-SNP subset) is that at least a bootstrap sub-sampling process is carried out, to make sure that observations found within the comparison of interest are at least well replicated when comparing smaller sub-groups of the cases and controls. More information on that can be found in the poster I presented at Queenstown Research Week last year, and detail about the implementation in my draft preprint (which will probably never be properly published, given how old the research is).


Reads mapped to exonic, intronic and intergenic regions

As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.


How can I assign a genomic region into a window using R?

I'm not sure how you did your calculation but for me it looks like your region is 3499713 bp in length.

site <- GRanges("chr1:6199548-9699260")

seqlength <- c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)

bin_length <- width(site)

## check length of region
print(bin_length)
##3499713

bins <- tileGenome(seqlength, tilewidth=bin_length, cut.last.tile.in.chrom=TRUE)

## extract overlapping bins
flank.bins <- subsetByOverlaps(bins,site)

## find new coordinates excluding region of interest
new.flank.bins <- setdiff(reduce(flank.bins),site)

## combine new bins with your region of interest
new.bins <- c(new.flank.bins,site)

## remove old bins
bins <- subsetByOverlaps(bins,new.bins,invert=T)

## add new bins
bins <- sort(c(bins,new.bins))

## check that it makes sense
##table(width(bins[seqnames(bins) == 'chr1']))

write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names= FALSE, row.names = TRUE)


Old versions of the reference genome and dbSNP

The UCSC Genome Browser downloads page has old human genome assembly versions way back to the first human genome ever published (by UCSC in 2000).
http://hgdownload.soe.ucsc.edu/downloads.html#human

We also have a variety of archived dbSNP data up to dbSNP141 on hg38, dbSNP 131 on hg37, dbSNP126 on hg18:
http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp141.txt.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/snp138.txt.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg18/database/snp126.txt.gz

They're in a custom table format, which can be turned into VCF manually.

You can also get all of this data from Table Browser and add filters, etc.:
http://genome.ucsc.edu/cgi-bin/hgTables


How do I efficiently subset a very large line-based file?

Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve.

Since I don’t know how how to do this in sed, and it’s not trivial in awk either, here’s a Perl script:

#!/usr/bin/env perl

use strict;
use warnings;

my $file = $ARGV[0];
my $lines_file = $ARGV[1];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;

open my $fh, '<', $file or die "Cannot read file $file";
my $line = 1;
my $next_line = 0;
while (<$fh>) {
    last if $next_line == scalar @lines;
    if ($line++ == $lines[$next_line]) {
        $next_line++;
        print;
    }
}
close $fh;

I’ve implemented a similar function in C++ for an R package, that's only slightly longer than the Perl script. It is ~3 times faster than the Perl script on my test file.


How to correct alpha, and not p-values themselves, for visualization purposes

The only way to get the alpha levels is to determine what they will be with p.adjust(), since they will depend on the distribution of your unadjusted p values. The general steps you should follow will be:

Add a column of adjusted p-values to your dataframe (mydata$padj = p.adjust(mydata, method="BH"), which is the same as FDR and saves a character).
Use which and max to determine your two alpha threshold (e.g., max(mydata$pvalue[mydata$padj < 0.05])

Then you can adjust your plots however you like (presumably with some horizontal lines at the various alphas). Whether you take the smallest non-significant value or the largest significant value is up to you, just describe what "dots on the line" represent.


Is there a way to tell which chromosome a gene is on, by looking at the "Chromosome/scaffold name"

What you are looking at is a patch, not an actual chromosome location. HG-2128 is an issue ID, something that the Genome Reference Consortium uses to track issue with references. Information regarding that specific patch is here. The issue is described as a "Possible misassembly or indel variation in GRCh38 within AL049860.8."

During the alignment step before you received your data, someone used a top level assembly version of the reference genome. This includes extra, accessioned scaffolds tacked on to the end of the file that contain the fix. These scaffolds represent what will be found in the next major release of the genome. The fixes aren't applied directly to the reference as soon as they are found due to the inevitable disruption of coordinates. So reads mapped to this FIX patch rather than the primary original scaffold in the primary assembly, on which the RNU6-280P gene resides.

The coordinates are slightly different. In the current GRCh38.p12
 assembly the RNU6-280P gene is located at chr6:67546651-67546754 (source). You can see that the coordinates are very close to those indicated in the table you provided (Gene Start, Gene End) but are indeed different.

Here is some more information on patches.


Your safest bet is to just redownload the annotation. I've not seen official documentation of the index used by annovar, but apparently it's just a text file with a single line header and chromosome-bin coordinates. Over on SEQanswers, there's a thread discussing this issue with a single perl script that allegedly recreates the index.


TPM or rlog(CPM) for comparing expression?

For comparing the counts of different samples from DESeq2, Michael Love recommends using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:

## [assuming "dds <- DESeq(dds)" has been run]
dds.counts <- assay(vst(dds, blind=FALSE));
## [This makes the scaled counts resemble more closely the actual read counts]
dds.quantile99 <- quantile(dds.counts[dds.counts > min(dds.counts)], 0.99);
dds.counts <- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);

The matrix dds.counts should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:

gene.VST <- sort(dds.counts["geneOfInterest",]);
plot(gene.VST);


At what processing step should library strandedness type be taken into account?

I usually take into account library strandness during the mapping to the reference step (but I never worked with bigwig files myself).

Assuming RNA-seq data from Illumina, you can use Hisat2 for alignment to the reference genome, and you have the option --rna-strandness to specify the strand-specificity information (default unstranded). This is particularly important when you wish to have XS strand attribute in your bam file for easier compatibility with downstream software (especially important in my experience if you are planning to use some utility from Cufflinks for your step 3).


Find a map/correspondence between two versions of a genome

A relatively convenient method would be to create a .chain file, as used by tools like liftOver and crossmap (and probably remap from NCBI). For your use-case, you can simply use the DoSameSpeciesLiftOver.pl command as described in the link. The resulting chain files are pretty self-explanatory and can be used to map BED files of regions of interest.


Gene expression signature (AUC) from gene expression data with 48 +48 samples

I have a method for doing this that uses a MCMC algorithm originally designed for population genetics (via the structure program) to generate a statistic between 0 and 1 indicating genetic risk. This needs the use of reference populations for cases and controls that the query individuals are added to for classification. I use T1D as an example. My methods paper is not peer-reviewed, and probably doesn't go into sufficient detail for what you are looking for, but at least indicates one way this can be done:

https://www.biorxiv.org/content/early/2017/01/31/104497

Based on the work I've done with this, using 48 cases and 48 controls with a very specifically defined type of cancer should be sufficient to create a generalised signature from human populations. If a generic cancer profile is needed, the populations are not representative of a general human population, or there are some cases in the control group, then things get trickier and require more individuals. On the other hand, 1500 per group is probably overkill.


DMPfold: Understand the run of DMPfold "run_dmpfold.sh"?

From the docs:

  The final model is final_1.pdb and other structures may or may not be generated as final_2.pdb to final_5.pdb if they are significantly different. Many other files are generated totalling around 100 MB - these should be deleted to save disk space if you are running DMPfold on many sequences.

These other files include intermediate CNS files (sequence_sub_embed_1.pdb etc.), output files before clustering (sequence_1.pdb etc.) and saved ensembles of previous iterations (ensemble.1.pdb etc.).

Let the run finish or you won't get the final output files. If you want to use all the generated structures, use the ensemble.n.pdb file from the last iteration (i.e. the highest number). You can also use the other ensemble files if you want to include models from previous iterations.

I would strongly recommend using default "3 50" parameters until you get things working. Don't move files during the run. You really don't get much benefit of going beyond 50 models per iteration, trust me. We tried it relatively extensively!


Can I eject the flow cell MinION from USB port when doing basecalling?

You can, you will lose the transfer from the MinION. It would probably be better to just stop the run. The computer will continue to basecall untill it is done. All the raw reads are stored locally on the machine so it should be fine.


WGCNA module preservation error code with plotting summary statistics

My first guess is that there is a typo in your code, and Acutal.IGP should be Actual.IGP. I checked my tutorials and the word Acutal is in none of them, so I would assume you mistyped the variable name.


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


Download program for ΔG for TM insertion

The application is stored in github (ElofssonLab/dgpred) with a link to the website.

To find it I searched for: dgpred download github


Comparing two files with overlapping regions and get their associated information


How to visualize called narrowPeak files in UCSC Genome browser or IGV?

I think you might have changed the separator (or at least have some kind of inconsistency from the required format) for your file. Note that peak output files from MACS2 are variants of BED files.

It seems you need to have a tab separator for this type of file.

I copied your example file and then ran awk -v OFS="\t" '$1=$1' your_peaks.narrowPeaks > fixed_your_peaks.narrowPeaks on it. This produced a result that I was able to see in both the UCSC genome browser and IGV.


How to make Venn diagram for the Macs peak calling output of two data sets?

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


What is the difference between fixed effects and random effects in the context of Linear-mixed models?

I do not know much about statistics but I will try my best to explain.

First, random effects are defined as the factors (categories) in the population that we are not aware of (not observed), so we are randomly sampling levels of those factors when we sample the population.

Practically speaking, random effects can be found when there are hierarchical structures in the data. For example, we want to compare the relationship between students' SAT score and their college admission rates. You might just want to run a simple linear regression. But there might be some structure in the data that is hidden. For example, the students can be grouped into classes, and then schools. You can also group students by the social-economic status of their parents.

The idea is that students belonging to different categories may have different college admission rates even if they have the same SAT score. Therefore, it is important to consider some of these "random effects" when building your model.


GATK 4.2 using IndexFeatureFile with multiple files

You could run loop over the files in your shell, running a separate command for each file. For example, in Bash (or another Bourne shell):
for file in *.vcf.gz; do ./gatk IndexFeatureFile -I ${file} -O ${file}.tbi; done


GSEA - Gene Set Enrichment Analysis - in R studio

I have answered my own question however, it is specific to windows users.. If you get the same error from the cmd as I did, you must manually put in environment variables related to Java. Do as the user Mohamed did in this SO post: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for

Following this, you must avoid using paths such as C:/absolute_path and instead use /absolute_path, if you do the former, you will get the exception: java.net.MalformedURLException: unknown protocol: c

In the end, after manually entering the environment variables and running the following command: java -Xmx2G -cp /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file /Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx /Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out /Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7 it has worked.

Where the command can be executed directly from the command line or in R using the system function.

Note: A '#' separated by a space must be present at the beginning of your .rnk file for the program to work. Therefore, .rnk files must have the following format:

# Genes  Coeffecients
TSPAN6    0.2345676
etc. for every gene.


Visualising gene expression across cell type and conditions in one plot, in Single Cell Sequencing data

Try using SplitDotPlotGG, this will allow you to color the dots by some additional variable (eg condition).

In your case you would want to do something like this:

library(Seurat)

SplitDotPlotGG(object = object,
               genes.plot = c("Gene1", "Gene2"),
               plot.legend = TRUE,
               group.by = "cell_types",
               grouping.var = "tissue_types")

This would give a result similar to below (an example from https://satijalab.org/seurat/immune_alignment.html).


Basically and briefly,

Normalisation = transforming a data set to conform to a normal distribution. Converting to a frequency is NEVER a normalised transformation
Standardisation = ensuring each point of the data set is comparable to the others, either directly or relatively depending on the calculation

The information you have received is incorrect, moreover you only convert to a frequency for a final assessment of the data and in my experience NEVER as an intermediate step. Normally data point between 0 and 1 are bad news and often require a transformation to get them away from fractions.
I don't understand the behaviour of the variance involved in metabolism and particularly lipids, however you should refer this issue to "cross-validation" Stackexchange. I did have some understanding of enzyme kinetics but in a different area.
Cross-Validation Stackexchance statisticians will understand less of the biology than I do (which is a bit worrying of course), so you will need to give them quite detailed information on the data set. They will almost certain support my concern of the advice you are given and provide a better method of standardisation. One method is to take the difference in mean and divide it by the standard deviation, you can then normalise the standardised data. They will want to know the goals of your investigation to understand who to process the data further.
What I can say is variation between runs is often an issue (in your apparatus perhaps not, I dunno) and standardising against a common control is one approach to overcome this, e.g. the control same common to all runs would be the mean about which all other data is centred.
What I suspect is that lipids are less popular than other fields, so standard protocol may not be observed.


Calculating most abundant transcript from RNA-Seq data

Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here:

salmon (as you mentioned) to quanitfy. Run it with --numGibbsSamples 50 (or higher if your computer has the hardware. Also run with -d for the next tool
terminus is a tool that also originates from the COMBINE-lab, It's relatively experimental, but it applies well to your use case. It will take the salmon output and assign counts to transcript groups. If a set of transcripts can be clearly quantified above a certain confidence, that transcript will stay unique. If it is confounded with other transcripts, those will be placed in a transcript group and their counts will be summarised. More info here

You are right in that TPM, FPKM etc. are not appropriate measures to compare across samples. TMM, or DESeq2's median ratio method (what you get with counts(dds, normalized=TRUE) are a step in the right direction, but are still inappropriate across samples. I'd suggest this often posted blog post in the topic.
There are no perfect solutions here since RNA-Seq is inherently a relative quanitifcation method so this will always be an issue. IMO, use DESeq2::varianceStabilizingTransformation(), DESeq2::rlog() or limma::voom() after terminus and rank your transcripts that way.
Just keep the inherent transcript incertainties in mind.


Programmatically retrieve Accession List from NCBI

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

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

to download the table in CSV format.