Calculating Isoelectric point from a multifasta file

You can use EMBOSS Pepstats for this. It takes a multi-fasta as input and produces a table that contains various statistics related to the protein, including isoelectric point.


How to obtain file offsets corresponding to range query in block-gzipped file

So if the purpose of obtaining file offsets is to access subportions of remote files via protocols that samtools supports, as mentioned in this comment, then you can make your life easier and use samtools/tabix directly on the cli, prefixing your datafile with the corresponding scheme.

One side effect of using the CLI however, is that the blocks retrieved will be decompressed (e.g. samtools view). You'd have to re-compress them if what you're trying to do is create a smaller "bam" file.

Keep in mind that if you retrieve all zipped chunks which cover your range of interest, then you might be extracting more data than you asked for. The first and last chunk may contain records outside your region, just because their neighboring records happened to be compressed in the same chunk. You'd have to be careful to filter those downstream.

Specification for the tabix format

If you absolutely need a custom solution, the tabix data structure format, authored by Heng Li, is on the samtools website:

https://samtools.github.io/hts-specs/tabix.pdf

A more in depth discussion of how the data structure works can be found here:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3042176/

The tree-shaped binning index is a bit hard to understand from the text alone, but the tabix paper mentions that their approach is almost identical to the binning algorithm in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC186604/pdf/X6.pdf (page 1003). That paper has a nice illustration which makes things clearer.

Operation

Note: This is from my understanding of how the indices interoperate.

You don't need the data file to retrieve chunk offsets, you can do the work with just the tabix file. But it's complicated, and probably error prone if you're implementing this yourself.

Let's assume the query interval is Chr01:100-Chr02:2000.

Each sequence (Chr01, Chr02, etc) has its own pair of indices (binning index + linear index) in the tabix file. So you would break the query into two queries: Chr01:100-999999999999 and Chr02:0-2000, and concatenate the two results.

Within one sequence, to find all the blocks for open interval [X,Y), you would use the binning index to retrieve all the bins overlapping this interval (There is C code to do this in the data structure pdf). Then, you need to go through all the possible bins selected, and filter chunks. Not all chunks in a given bin will be overlapping with X, Y.

To filter chunks within a bin, you need the help of the linear index. There is one entry in the linear index per 16kbp of sequence positions. Each entry of the linear index contains the virtual offset of the leftmost record overlapping that 16kb window. I don't know what the entry contains if there is no record overlapping that 16kb window, presumably 0.

You would find the linear entry in which position X falls (X/16Kb), take that offset, call it O_x. You find the linear entry corresponding to position Y, and take that offset, call it O_y. When you go through the selected bins, filter out any chunk whose offset is < O_x, or whose offset >= O_y.

Assuming your datafile is sorted, the chunk range you're interested in is contained between the minimum offset and the maximum offset selected.

The virtual offsets are 64bit unsigned integers. The higher 48 bits of the virtual chunk offset corresponds to the file offset of the datafile (the other 16 you don't need. They are the offset of the record within that decompressed chunk).

Software

The software I have found typically mimics closely the hts_lib API. And hts_lib doesn't export a function which loads an index without also requiring a datafile. As far as I know.


Convert rs ID of one hg build to rs IDs of another build

You can assume that the overwhelming majority of rsIDs are the same between GRCh37 and GRCh38 (they're semi-stable IDs). There are, however, a number of rsIDs that are present only in GRCh37, which you can find here. Note that the format of this file is a bit strange, it's chromosome|position|ID|weight, where position is sometimes empty, weight is typically 1, and ID is the rsID without the rs prefix.

Also, some rsIDs have a different strand. This is less likely to be an issue for you, but it's there in case you need it.


t-sne code and labelling issue for rna seq data

Although a bit late, here is the solution for this:

1) Create color/column levels:

levels = unlist(lapply(colnames(trn[,1:16]), function(x)
    {

        if(x == "H1") return("H1")
        else if(x == "H2") return("H2")

}))

2)

plot(tsne$Y, col=factor(levels))

This should color columns H1 and H2 differently

Edit:

This is a sample complete code:

library(Rtsne)
file1<- read.csv('filename.csv',header = T, row.names=1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=5, theta=0.5, dims=3)

collevels = unlist(lapply(colnames(trn), function(x)
{
  if(grepl("H1",x))
    return("blue")
  else if(grepl("H2",x))
    return("red")

  # etc ... for all the columns you want to color
}
))

colorColumns = factor(collevels)
require(rgl)
plot3d(tsne$Y[,1], tsne$Y[,2], tsne$Y[,3], col=colorColumns)


Note that in the article published that figure is different from the preprint!
See the figure:

As they aren't substantially different in display other than showing the two set of predictions that are comparable between CAFA1 and CAFA2 I explain them below:

The colors are the margins of $F_{max}$ they refer to the equation found in the section "Comparison between CAFA1 and CAFA2 methods" under methods:

$$\delta(m_2, m_1) = \frac{1}{B}\sum_{b=1}^{B}F_{max}^{b}(m_2)-\frac{1}{B}\sum_{b=1}^{B}F_{max}^{b}(m_1)$$

Where B is the 10000 bootstraps they did and $m_1$ and $m_2$ are the CAFA1 and CAFA2 methods.

Where $F_{max}$ is explained in the article (and in the preprint) as:

  A perfect predictor would be characterized with $F_{max} = 1 $

So, the higher the (positive) difference is, the better precision and recall (See wikipedia) are for those methods. Note that the scale is from -0.1 to 0.1 not from -1 to 1 so the differences are quite small (One could argue that the log10 would be better to see those small differences)

One can appreciate that in the bar plots at the right, where the differences between each dataset are small in the B section, so the color is white, and in the C panels they have brighter color.

The numbers in the boxes are the percentage that a CAFA2 method did better than a method in CAFA1 for another dataset in the "intersection of the two GO snapshots" as noted in the methods section. (But it is also not clear to me)


Getting a stretch of genomic ranges from a dataframe/granges object based on metadata column

The GenomicRanges package in bioconductor will do this.  It lets you do 'set' operations on genomic ranges, such as concatenating them, finding overlaps, etc.  Hopefully something in the following code will be inspiring, even if I haven't solved the exact problem for you here:

library(GenomicRanges)

# make a GenomicRanges set corresponding to the original table
test.gr = makeGRangesFromDataFrame(df = test, start.field = "start", end.field = "end", seqnames.field = "chr")

# make separate GenomicRanges sets for each of the two partitions of the genome
test_A.gr = makeGRangesFromDataFrame(df = test[test$comp=="A",], start.field = "start", end.field = "end", seqnames.field = "chr")
test_B.gr = makeGRangesFromDataFrame(df = test[test$comp=="B",], start.field = "start", end.field = "end", seqnames.field = "chr")

# apply the reduce() function to concatenate segments which are adjoining each other into single segments
reduced_test_A.gr = reduce(test_A.gr)
reduced_test_B.gr = reduce(test_B.gr)

# group PC1 values by segment classification
test[queryHits(findOverlaps(test.gr, reduced_test_A.gr)), "PC1"]
test[queryHits(findOverlaps(test.gr, reduced_test_B.gr)), "PC1"]


Why does `pratchet()` in `phangorn` not find the most parsimonious trees?

Perhaps the tree which you have calculated is not being caught in the search heuristic? Have you tried (as a trial) setting the tree which you calculated as the starting topology, just to check that it is indeed being calculated on by pratchet()?

Alternatively, you could try playing around with the search heuristic settings, like varying from ratchet to NNI, or you could also try using an external program, like PAUP*.


Merge / Reconciliate several de novo transcriptome assemblies with different kmers

You may combine all your transcriptomes into a single file and then apply a clustering method to group very similar transcripts into a single one. For that purpose, you may try CD-HIT-EST or MMseqs2. For each identity threshold you are going to test, you may assess the final quality with BUSCO or by blasting against reference sequences.


Suffix arrays vs Suffix tree is a problem quite often discussed in bioinformatics classes, but it's not that much used in practice. Most of people use already optimized mappers bowtie2 or bwa-mem.
I would advise you against your intention of mapping the reads using your own implementation of anything. You don't need to reinvent the wheel, just use an existing solution.
By the way, the Ns in the reference are either hard-masked nucleotides (that would be repeats of some sorts) or (more often) simply unknown nucleotides. So there is no way to map your reads there anyway (and the mentioned mappers will handle the Ns for you).


Can structural - equation modeling be used in Bioinformatics?

I have came across them in gene regulatory network analysis before. see here and here. However I haven't seen them in other contexts, but that's not to say they don't exist.


Downloading multiple SRA files from several SRA accession IDs does not work

First of all, the brace expansion you tried to use isn't valid. While bash can expand {1..3} to 1 2 3, it has no way of knowing how to expand {SRR1002678..SRR1184123}:

$ echo {1..3}
1 2 3
$ echo {SRR1002678..SRR1184123}
{SRR1002678..SRR1184123}

Perhaps you meant SRR{1002678..1184123}? But there's no need for this, you already have the accessions in a file and can read that. I see that the FTP server you are accessing has organized the data like this:

sra/sra-instant/reads/ByRun/sra/SRR/SRRNNN/SRRNNNXXXX/SRRNNNXXXX.sra

Where NNN are the 1st three digits of the accession. So, for example, the file for accession SRR1002678 is at:

ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR118/SRR1181152/SRR1181152.sra

With this in mind, you can write a little script that extracts the URL from each accession:

while read i; do
    n=$(echo "$i" | sed -E 's/(SRR...).*/\1/');
    wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/$n/$i/$i.sra
done < SRR.txt

The only trick is the sed command which simply looks for SRR and the three next characters and removes everything else on the line (it "captures" the found pattern and then replaces the line with what was captured (\1)).


I could install HiC-Pro using conda when creating a separate environment (which I always do for side-projects), and then pull the conda build of HiC-Pro:
conda create --name HiCPro python=2.7
conda activate HiCPro
conda install -c davebx hicpro

Not sure what conda.sh is that you describe, sometimes tricky to get random code from GitHub running. Give it a go with a standard conda install as I described.


SNP located within a promoter region (pig)

As far as I'm aware, Illumina provide CSV annotation files for all their sequencing chips, which can be used when they can't be found in Bioconductor. You can find annotation information for the PorcineSNP60 here, in particular the Manifest file (CSV format). The format is Illumina's weird "we say it's a CSV because there are commas in it" format, so if there are only a few locations of interest it would be easier to use grep than try to work out the right way to load it into R (or similar).

For completeness, you want the lines between [Assay] and [Controls]. Here's an example way to pre-process the "CSV" file into a more standard format that just has one data table in it:

$ grep -n '^\[' PorcineSNP60v2_15031945_C1.csv
2:[Heading]
7:[Assay]
61574:[Controls]
### See that the table starts at line 8 (after [Assay]),
###   and finishes at line 61753
$ tail -n +8 PorcineSNP60v2_15031945_C1.csv | head -n $((61574-8)) > cleaned_PorcineSNP60v2_15031945_C1.csv

But for a few markers and grep, it doesn't matter if you're working with the clean or dirty CSV file:

$ grep -e '^IlmnID' -e '^MARC0073381' -e '^ALGA0066960' PorcineSNP60v2_15031945_C1.csv
IlmnID,Name,IlmnStrand,SNP,AddressA_ID,AlleleA_ProbeSeq,AddressB_ID,AlleleB_ProbeSeq,GenomeBuild,Chr,MapInfo,Ploidy,Species,Source,SourceVersion,SourceStrand,SourceSeq,TopGenomicSeq,BeadSetID,Exp_Clusters,Intensity_Only
ALGA0066960-2_B_F_2199219285,ALGA0066960,BOT,[T/G],0054749354,TGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA,,,10.2,12,55218591,diploid,Sus scrofa,PorcineSNP60v2,2,BOT,GGCTCACCTCTGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA[T/G]CCAGGCTTCTCTCTCCTACTCCTAGGGCCCCTGATGGTTCCTGCATCCTGACCAATAGTG,CACTATTGGTCAGGATGCAGGAACCATCAGGGGCCCTAGGAGTAGGAGAGAGAAGCCTGG[A/C]TAGACCTGTATCTTGCGAGTGTTGGGCTTGAGCTGAGCAGACCTGTGGCAGAGGTGAGCC,670,3,0
MARC0073381-2_T_F_2199252059,MARC0073381,TOP,[A/G],0047761415,TGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA,,,10.2,7,61808556,diploid,Sus scrofa,PorcineSNP60v2,2,TOP,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,660,3,0

I can see from this that the first marker is Chromosome 12, position 55218591, and the second marker is Chromosome 7, position 61808556. The Illumina annotation doesn't include the nearest gene, so you'll need to hunt for those locations on a genome browser (or using mapping with Bioconductor).


Filtering bases based on phred qualities with pysam

The PileupRead object as a query_position attribute, which you can use for this:

for pileupcolumn in samfile.pileup() :
    for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
             continue
         if not pileupread.is_del and not pileupread.is_refskip:
             if pileupread.alignment.query_qualities[pileupread.query_position] < 10:
                 # Skip entries with base phred scores < 10
                 continue
             if pileupread.alignment.is_reverse: #negative
                   ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
             else:
                   ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

Note the 6th line, which implements a filter with a threshold of 10.


Normalization for two bulk RNA-Seq samples to enable reliable fold-change estimation between genes

You have only two samples?

You aren't going to be able to draw strong conclusions from that no matter what you do.  Clever statistics don't work without replicates.


How do I set a neural network to loop multiple times and average the resulting values?

The concern with your approach is called 'leakage' because you are parameterising the same data set that you are training. This can easily lead to overtraining is in ANN (artificial neural networks) and is not cool. Overtraining is where the ANN very tightly wraps around the training data set that it performs poorly on any other data set. In other words it is overtightened.
The way they get around it is to split the data set usually 60:20:20 into a training split, parameterisation (validation) split and test (testing) split. If the parameterisation split gets mixed with the training set - leakage happens and overtraining results. Normally the data is trained, and the trained model parameters are assessed/optimised in the parameter split. The trained data and optimised parameters are then tested in the test split.
As long as you are aware of the approach you are using, that I guess that is ok. If you separate your training and parameterisation (validation) data that would be cool.


How can I get or create a reference genome for Bacteria?

How can I get a reference sequence for other strains of E. coli?

You can use any sequence you want as reference for a variation analysis. It just depends on the question you want to answer.
So to keep it simple, your Question is:
What are the variant region from a fresh sequenced E. coli K12 sample?
To know the differences from your known strain E. coli K12. A reason for this could be that you work with the strain K12 through out several years. During your studies you changed the conditions for the bacteria, like adding antibiotics or change the growth media. So some of the E. coli would adapt to the antibiotic and become resistant to it or they adapt better to the new media. This happens because some nucleotide changed. With each new experiment you used the best adapted colonies from the previouse one.
So you sequence your bacteria and compare it with the E. coli K12 strain. In this case this is your reference. If you use another strain, like SE15 or SE11, as reference, your result would not give you any answer for your question.
If your used data and publications gives you the information which strain were used, then use the given strain as reference.


Stacked barplot for single cell analysis

Based on your example, you can use count for each phase of cell cycle for each cell types in each condition. Here, I'm using dplyr package to do that but most likely, you could have the same output using various methods:
DF <- rbind(dfhead, dftail)

library(dplyr)
DF_Count <- DF %>%group_by(orig.ident,Phase,old.ident) %>%
  count() %>%
  ungroup() %>%
  group_by(orig.ident,old.ident) %>%
  mutate(Freq = n/sum(n)*100)

# A tibble: 8 x 5
# Groups:   orig.ident, old.ident [5]
  orig.ident Phase old.ident      n  Freq
  <chr>      <fct> <fct>      <int> <dbl>
1 treated    G1    Fibroblast     1 100
2 treated    G1    Stellate       1 100
3 treated    G2M   Tumor          1  33.3
4 treated    S     Tumor          2  66.7
5 treated    S     Mast cell      1 100
6 untreated  G1    Acinar         1  16.7
7 untreated  G2M   Acinar         2  33.3
8 untreated  S     Acinar         3  50

As you can see, DF_Count has the frequency for each cell types of each phase of the cell cycle in function of the condition. We can use DF_Count to get the following plot. Using facet_wrap, you can create two panels based on the condition column and thus represent the cell cycle of each cell types in function of the treatment condition:
ggplot(DF_Count, aes(x = old.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~orig.ident)

Does it look what you are expecting ?
NB: DF here is the addition of your head and tail of your all.combined dataframe.


Why are some of the gene peptides returned by biomaRt missing an asterisk in the end?

If it's 3' incomplete that means the evidence used to create it was a fragment. Here's the evidence used to construct BRCA1-214 ENST00000477152.5, a 3' incomplete.

You can see that there's a full length cDNA from EMBL, AK307553.1, which was used to create this model.

The sequence was mapped against the genome sequence to create the transcript model. When searching for the translation, it was found that the open reading frame was started, but there was no stop codon. This suggests that the cDNA AK307553.1 is actually a fragment, that the mRNA was broken or cleaved before it was reverse-transcribed sequenced. We display it in Ensembl in the hope that this will lead someone to identify the full length transcript that it represents.

Many people choose not to work with these incomplete transcripts. You can do this by filtering by Gencode Basic in biomaRt.


You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file.

bedtools genomecov -d -ibam $inputfile > "${inputfile}.genomecov"

You can also provide a BED file if you just want to calculate in the target  region.


So far what i know is i have to perform a parametric test .

How did you determine this?

  But im not sure if I have to go for a

  Paired t-test or Unpaired t-test

If you can't tell the difference between when to do a paired t-test and an unpaired one, you need to stop and talk to someone who can teach you basic statistics.  Throwing statistical tests at numbers when you don't understand them is a recipe for disaster.


I remember having this problem in the past, it was encoding related. Rewriting the -s manually (the dash sign) should fix the issue, please see this.

  Problem is a copy & past from the BLAST documentation, all dashes
  ('–') need to be replaced by hand at the command line.


Early MinION sequencing runs had forward and reverse DNA templates joined together by a hairpin adapter, so that the sequencer would read both strands from the same template. The consensus sequence that was generated by combining the base calls of opposite strands is referred to as a 2D read.

Due to a dispute with Pacific Biosciences, and for other reasons related to chemistry and basecalling, the equivalent mode in current sequencing runs is called 1D². There is no physical ligation of forward and reverse templates, but the kit and flow cell chemistry are modified to encourage the reverse strand to follow the forward strand very soon afterwards through the pore, with the joining and consensus generation carried out in software after the sequencing is done. In the output fast5 files, these 1D² reads can also be referred to as 2D.


Parsing gene expression value for a list of genes from GEO

GEO provides a GUI for this (just click on the "Analyze with GEO2 button"), you can even check/download the R/Bioconductor code used for the analyses.
Basically you will create groups of samples to compare/contrast and R/Bioconductor will run behind the scenes. Here is a nice tutorial.


The tutorial you linked to does not mention ERCC or mt genes.

Before running grepl or any other search, you should first check your inputs. It looks like you have WormBase gene IDs if all your genes look like WBGene00023193. Therefore, your question is actually how to convert WormBase gene IDs to gene symbols.

Here is one potential solution: https://github.com/WormBase/website/issues/4297

After you convert, you have to make sure that you have ERCC spike-ins in your gene table (probably not) and you mitochondrial genes start with "mt-".


Tools for modelling and visualizing growth of cells

If the scenario is simpler, and you can't find anything by searching, I would consider just writing your own.

The modelling can be done in any language but if compute time is an issue - for example you want to run a simulation with many cells or many simulations - then I would use a language that can support that such as Julia or NumPy in Python. In fact Julia is a nice language for this type of thing as you end up with readable code that runs fast.

For visualisation in Julia specifically you could use the animation functionality of Plots.jl, which is basic but there, or could try Makie.jl. See the Molly.jl documentation for some examples of visualising particle-type processes using Plots.jl.


Publicly available genome sequence database for viruses?

There area few different influenza virus database resources:

The Influenza Research Database (IRD) (a.k.a FluDB - based upon URL)

A NIAID Bioinformatics Resource Center or BRC which highly curates the data brought in and integrates it with numerous other relevant data types

The NCBI Influenza Virus Resource

A sub-project of the NCBI with data curated over and above the GenBank data that is part of the NCBI

The GISAID EpiFlu Database

A database of sequences from the Global Initiative on Sharing All Influenza Data. Has unique data from many countries but requires user agree to a data sharing policy.

The OpenFluDB

Former GISAID database that contains some sequence data that GenBank does not have.

For those who also may be interested in other virus databases, there are:

Virus Pathogen Resource (VIPR)

A companion portal to the IRD, which hosts curated and integrated data for most other NIAID A-C virus pathogens including (but not limited to) Ebola, Zika, Dengue, Enterovirus, and Hepatitis C

LANL HIV database

Los Alamos National Laboratory HIV database with HIV data and many useful tools for all virus bioinformatics

PaVE: Papilloma virus genome database (from quintik comment)

NIAID developed and maintained Papilloma virus bioinformatics portal

Disclaimer: I used to work for the IRD / VIPR and currently work for NIAID.


PASA pipeline: compare experimental transcripts to the reference annotation

Answer from the authors:

  PASA ends up using gff3 format instead of sam format for uploading the
  alignments.  If somebody has the gmap gff3 output, then can be done to
  just upload that directly using the custom alignment importer, but not
  with sam. They suggest just having PASA rerun GMAP as part of its
  regular routine.


Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the read.has_tag('XS') continue block.

I think that if you filter out is_secondary and is_supplementary you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly.

for read in bam.fetch()
    if read.is_secondary=true
        continue
    if read.has_tag('XS') <- This block stores XS and then moves on to the next read
        XS = read.get_tag('XS')
        continue
    if read.has_tag('AS')
        AS = read.get_tag('AS')
    if AS < XS  <- This block compares the AS of the current read to the XS of the previous read that had one.
        read.write(read)
bam.close()


Merge hundreds of small BAM files into a single BAM file

samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO quickly becomes a bottleneck anyway.


How to calculate p-values for fold changes?

Your null hypothesis would be that the fold-changes are 0, so you can either do the T-test accordingly or simply do away with the fold changes and perform the T-test between the raw values in the group (this is preferable to performing a T-test of one group vs. 0, since it allows you to assess the expected variability around 0). Note, however, that 3 samples is pretty much the bare minimum needed for any statistics, so your power will be terrible and your estimation of variance will likely be pretty inaccurate. Consequently, take any p-values with an appropriate grain of salt.

To my knowledge there aren't any good alternatives for your situation without having a good expected background distribution.


Questions regarding the detection of variants (SNPs and Indels) in mammal

This is close to what many people do. May I add:

between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this.
filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold.
For that purpose, some time ago I made this and this.


If the tag r is incorrect, the easiest way is to download and install PyMOL: https://pymol.org/2/
Drap and drop your pdb file and type remove not chain A.
Then save molecule via the menu as pdb or type save some_file_name.pdb (will be in your home root folder).


Yes, you can use limma for this mixed model approach. Like you suggest, the random effect (persons) can be put in duplicateCorrelation().

Here is a similar example with RNAseq data, on bioconductor support site.


tool to visualize a collection of phylogenetic trees

I recommend using densitree: https://www.cs.auckland.ac.nz/~remco/DensiTree/

And rwty: https://github.com/danlwarren/RWTY

Densitree plots the superimposed set of trees in a heatmap-like format:

And rwty, a very powerful R package that allows visualisation of treespace in an NMDS format (among other methods):


Why this ballgown error (first column of pData does not match the names of the folders) comes out

Well, it looks like if you're using Windows for managing files - that's not a good idea - it sorts names of your ballgown folders in alphabetical order while you need it differently.

That is why the one should use Linux - when you're ls -ltr your ballgown folder in Linux it gives you the "right" sequence of folder names to be put inside a ballgown object using pData.

This is what helped me to get rid of this error:

In Linux enter your ballgown folder and do ls -ltr | awk '{print$9}' and you will have the right sequence that you have to put inside you .csv file using Notepad++, indicating what's what. Below's my example compared to the above mentioned initial example:

id,state
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube

I don't have my error anymore and after

phenodata <- read.csv('my_phenodata_file.csv')
bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata)

my results are nicely wrapped

Sat May 02 22:11:57 2020
Sat May 02 22:11:57 2020: Reading linking tables
Sat May 02 22:11:59 2020: Reading intron data files
Sat May 02 22:12:04 2020: Merging intron data
Sat May 02 22:12:06 2020: Reading exon data files
Sat May 02 22:12:14 2020: Merging exon data
Sat May 02 22:12:16 2020: Reading transcript data files
Sat May 02 22:12:20 2020: Merging transcript data
Wrapping up the results
Sat May 02 22:12:22 2020

The problem is explained here as well

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

Hope this will help anyone


How do i map subcellular localizations to lower the number of classes to predict

One way to do this would be to plot the first 10 GO categories of the the CC sub ontology (CC from cellular component). Or you can just visit the official webpage to see the list of children terms of them.

However there aren't only 4 or 5 categories, so you could arbitrary pick some of those.


Are there any databases for gene co-expression or expression pattern clustering?

How about:

http://www.informatics.jax.org/expression.shtml
http://www.emouseatlas.org/emage/home.php

Other potential databases of interest:

https://biokeanos.com/search?q=co-expression+%2Bmouse


How does DeepVariant construct RGB images from DNA sequences?

Actually the paper has made it clear how they did it. You just have to read the supplementary materials closer.

In their figure from the paper: A is Red, C is Green, G is blue, and T
is Yellow (G+R), but this is still unclear how they the 3xNxN image.

In RGB, each dimension is an NxN image. Since you have three dimensions, so it's 3xNxN. The red dimension was used to encode the nucleotide bases. The green dimension was used to encode quality scores. Finally, the blue dimension was used to encode the strand information.

Obviously, DNA is a string over an alphabet with the characters: {A,
T, C, G}.

This is easy. Lot's of ways. You could do a one-hot encoding or what DeepVariant used:
def get_base_color(base):
    base_to_color = {'A': 250, 'G': 180, 'T': 100, 'C': 30}
    return base_to_color.get(base, 0)


Find most similar genes to a set of genes of interest in time series RNA-seq data

You can use GSVA to summarize the information of your 8 genes and look for the time-serie change using limma or DESeq2.

If you have enough samples for each time point you could also do a weighted gene co-expression network analysis (WGCNA), and see if those clusters in each time point are co-expressed similarly (But it can be tricky and you will need to define a lower bound of the clusters to 8 or so). You could also divide your time series in two point at the beginning and at the end, and see if there are differences or not.

If you want to get a list of genes similar to your set of genes you can use GOSemSim or BioCor (Disclaimer I am a contributor to the first, and author of the second), which can calculate similarities between genes based on the gene ontology and on pathways respectively. For this approach you'll need to calculate the similarity between all the genes and then select those that are closer to your genes of interest.

I'm not aware of any database/website with time-series profiles of genes (I assume that if there would be any would be on humans or some simple model organism).


You need to follow the readme of this folder to  use the pipeline, there you will find instructions about how to align, mark duplicates, quantificate, and aggregate the results.

If you want to execute the WDL scripts you can follow this guide using Cromwell. Basically use: java -jar Cromwell.jar run myWorkflow.wdl --inputs myWorkflow_inputs.json


Instead of the "gene_id ..." bit, you want sprintf("gene_id \"%s\"; transcript_id \"%s:%s-%s\"", $4, $1, $2, $3).


Nextflow: Runs without error, produce few but not all output files into output channel

I suspect your HRC_1000G_check_bim and hrc_sites_file channels are queue channels when they should be value channels. Please again see how multiple channels work. Also, third-party scripts, like your Perl script, would be better placed in a directory called bin/ at the root of your project repository. Nextflow adds this directory to your $PATH, so scripts in this directory will be automatically accessible in your pipeline without needing to localize them using a separate channel. This way also avoids having to specify an absolute path to invoke them. Just make sure your script is executable: chmod +x your_script.pl
hrc_sites = file( params.hrc_sites_file )

process step5 {

    input:
    tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
    path hrc_sites

    """
    your_script.pl \\
        -b "${prefix}${chr_name}.step4.bim" \\
        -f "${prefix}${chr_name}.step4.frq" \\
        -r "${hrc_sites}" \\
        -h
    """
}


Download multiple fastq files using fastq-dump

A sample code is given in the salmon documentation as follows. Source
#!/bin/bash
mkdir data
cd data
for i in `seq 25 40`;
do
  mkdir DRR0161${i};
  cd DRR0161${i};
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_2.fastq.gz;
  cd ..;
done
cd ..

This could be modified as follows.
#!/bin/bash
mkdir data
cd data
for i in `seq 14 17`;
do
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_2.fastq.gz;
done
cd ..

You can save the code as a shell script and run it from the linux terminal. For example bash download_fastq.sh


NOTICE:
I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same.

This is something I have been meaning to get around to for a while, so thanks for the prompt.

I have created a python program called fast5seek to do what (I think) you're after.

As you have mentioned this is for educational purposes I have added a tonne of comments to the code too so I think you shouldn't have any issues following it.

The docs on the GitHub repo have all the info, but for those reading along at home

pip3 install fast5seek
fast5seek -i /path/to/fast5s -r in.fastq in.bam in.sam -o out.txt

What it does is read in <in.fastq|in.bam|in.sam> and extract the read id from each header. It then goes through all the fast5 files under /path/to/fast5s and checks whether their read id is in the set of read ids from <in.fastq|in.bam|in.sam>. If it is, the path to the file is written to it's own line in out.txt.

If no output (-o) is given, it will write the output to stdout.

So if you wanted to pipe these paths into another program, you could do something like

mkdir subset_dir/
fast5seek -i /path/to/fast5s/ -r in.fastq | xargs cp -t subset_dir/

The above example would copy the fast5 files that are found in your fastq/BAM/SAM to subset_dir/.


Can gene co-expression networks be used to help identify differentially expressed genes?

While you can use networks to find differentially expressed genes (see the WGCNA package, which does this) in my experience this ends up largely matching what you'd get using a traditional package with a looser threshold for significance. Given the time savings of traditional packages, there's rarely any gain to using networks (it won't hurt to try, just note that it'll take some time). If some very interesting genes are just above your significance threshold then change your threshold. You need to validate your findings in some way anyway, so your p-value threshold is partly just a way to protect you from wasting your time (but if background knowledge suggests that your time wouldn't be wasted...).


Which data is being used for violin plot?

Those two graphs don't look radically different.   The y axis is definitely not log scaled on the top one, and it looks like it is scaled in the second one (y.log = TRUE)


What are the de facto required fields in a SAM/BAM read group?

The sample tag (i.e. SM) was a mandatory tag in the initial SAM spec (see the .pages file; you need a mac to open it). When transitioned to Latex, this requirement was mysteriously dropped. Picard is conforming to the initial spec. Anyway, the sample tag is important to quite a few tools. I would encourage you to add it.


REMC SNP annotations for transcriptome prediction - epiXcan

After much research I found out how to obtain the annotations.

You must go to the REMC website: https://egg2.wustl.edu/roadmap/web_portal/

Go to the chromatin-state-learning tab, and choose the n-chromatin state model desired - the core has 15 states, then 18, then 25 (imputed states).

Following this, you must go to the download URL link and, you will find a lot of files, primarily with the following coding information: E0*.. This corresponds to a particular tissue/cell type and so you must choose the one relevant to your research question. The way they are correspond to each tissue/cell type can be found in the following paper:

https://www.nature.com/articles/s41467-019-11874-7#Sec22


What will be an appropriate mathematical distribution for SNP data?

Simulating genotypes with realistic correlation structures is indeed not so simple, and there's quite a few papers dedicated entirely to that (e.g. https://bmcgenet.biomedcentral.com/articles/10.1186/s12863-015-0173-4). Also, DEPICT (https://data.broadinstitute.org/mpg/depict/index.html) comes with a number of simulated GWASs to generate the nulls, so that's a simple way to obtain some ready-made fake data.

On the other hand, indepedent SNPs should behave as draws from $Binomial(2, MAF)$. Note that the standard QC procedure of filtering out genotypes that deviate from Hardy-Weinberg equilibrium is just a goodness-of-fit test against $Binomial(2, \hat{MAF})$. If you employ this filter, some LD-pruning, and still get different results, consider posting more details about the discrepancies you see - that would be indeed unexpected.


In your ATAC data frame, you have chr1   16103  16348. But in your merged data frame, you have chr1   16104  16348. Is it possible that you a is 1-based whereas your ATAC data frame is 0-based? This can happen in bioinformatics and can be annoying. This might be why you have a lot of NA in your merged data frame. You should subtract 1 from the start position of your a data frame while keeping the end position unchanged. Then try to merge again.

You should always be aware of your data format to know whether it is 1-based or 0-based.


What type of information does an .sf file contain?

You should provide more context, otherwise one can only guess what is it you are looking for! In case you are asking about Salmon's output, below is from their documentation:

  Salmon’s main output is its quantification file. This file is a
  plain-text, tab-separated file with a single header line (which names
  all of the columns). This file is named quant.sf and appears at the
  top-level of Salmon’s output directory. The columns appear in the
  following order:

Name   Length   EffectiveLength   TPM   NumReads


Intron retention events for a condition sample with 3 replicates

All introns will have some level of retention and I guess no intron will be 100% retained. In the absence of a comparator, you will need to work out how much intron retention an intron needs to be called retained. This is a human judgement, there is no automated way of deciding this.

What you need to do this is some sort of PSI calculation (percent-spliced-in). MISO will provide this, but it doesn't do replicates, so I guess you'll need to run it on each sample and then average the PSIs and determine if it surpasses some sort of threshold you have set.


If you want to stick to grep, use a scripting language such as Perl to generate the regex programmatically. For example:
perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' ATCGCTATCG

Prints:
[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]

You can use it in grep like so:
grep '[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]' <<< ATNGCNATCG

Prints:
ATNGCNATCG

If that works for you, you could make it into a little bash function that also runs the grep. Add these lines to your ~/.bashrc:
grepN(){
    seq="$1"
    file="$2"
    pattern=$(perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' "$seq")
    grep "$pattern" "$file"
}

You can now run
grepN ATCGCTATCG my.fastq

Of course, this is not a good idea since the sequence might be in different lines, but that's what you were doing originally.


[disclaimer: not a PacBio person, I'm just vaguely aware of the technology]
Movie files are the PacBio raw data format, representing observed signal intensities in each sequencing well over the course of the sequencing run. The movie files can be used to re-call sequences with updated statistical models, which can improve accuracy, or test out hypotheses about how epigenetics alters the rate of nucleotide incorporation (i.e. methylation).
Illumina has a similar static version of this which is an image (or multiple images, depending on chemistry) of the flow cell after each synthesis cycle. Illumina's raw data format consumes tens of terabytes per run, and is of minimal value after base calling, so is almost always discarded after base calling.
Anything uploaded to NCBI SRA must contain FASTQ information. While the original uploaded files can be kept (which may include time/signal information), the files are converted into a FASTQ-like format for NCBI's internal representation. As such, there's a chance that movie files could be uploaded to SRA, but those files would also need to include the called FASTQ information (e.g. in a HDF5 container). In contrast to this, ENA does allow uploading raw data files without called information.


RefSeq transcript naming - transcript vs gene

No, this is the same sequence. A single RefSeq accession will correspond to one sequence. However, nobody said that sequence needs to be unique in the genome. The specific sequence you are talking about, NM_001085545.1, is a mouse transcript which can be aligned to many different regions of the mouse's chromosome X.

This is your sequence:

>NM_001085545.1
agagaagtctaggcaggactggcagactacagcaatgcctcaaccaagaagcctctgcga
gagcttgttctaggaattttcatcaatggggcttttagtcagcagggtcttgagatgcag
ggattccagtctgttagagccacagccagaagccatagccggagccagctacattcctgg
cagtcgcaagcgaaaaagaaacagtttggaggagttggcaacaagttctaatgttcatgg
acctcaaaaccagggaatgtatccacatcaagttctcaactacatctactggaaaagggt
taagatctcatctaatgatgcttatcaaaacttatttttggacgggcatgatagcgacat
taaaatccgtgctctgggaaaaacatggtgtttacacaaagtatttttatgtcagtcagg
ctactttgctaacattctcaaaggtacttggagagaatcacaccatggtgttataaatct
gatcattaagaatgaggatattgatacccgatctctgcattttgtgtttggtgctttgta
cacggatgcggatttgtcaataacacctctggaagttcctcaagttttggcagcagcatg
cctgcttcgggtggatcaagtaattcagcagtgtgaaggaatcatgaaagaaactatcaa
caggaacactgtgtgctcctattatttggcagcagaaacctatagattaaaagctgtaaa
gacgagatgctttgaatggcttctttgcaatttgatggtacatccaagtgtggcacttta
caaggaagtagacttgaagttgatgtatcttctagcactgtcttctgacttactagtcat
gcaaaaggagattgatgtatataccacactaaaaatatggatgttcctttatcttaatcc
atgctggaacggaaccatgaaacagcttttacaacacgcaaacaactggctttccaccca
catggcatatgttgataacatcagttttcttgaaagtgaagaaggactaatatttcaacc
agtgtttaaaaagctgagatttcagcacattatctgtgacttgacttccacaactattct
tgaacaagatcgactaatacctatggcatggttgtcacccatttacaaacaacagtggtt
gactttgctgcgaacacaagaatatggggtaattggaccacaagttatcaatgaacaaga
acttgaagaatgcaccatgaggtgtggtacaatgatccccaaggatggaagatatacttg
gaagtggtcagttggacgacttggctttcctttacgtgtgacctttaccaggcagtgtgt
aattttaaggcaacggtgtcagaggtgtgatggttctgcttgccacaaccatatccgaaa
tgtcatattcagaataactttggtgtgttttgattccaacaaaagagtaactttcagaaa
gacaacaggttataaaatcctcacctttgaatataaggaggagcaaattgtaatgaaatt
ggatagtgatgttctaaccttccctatgtgtatattctgcaatttcctttttgtaaacct
aggaaatgcagaaaacaagtaatcttatcaattattagactgtaagtattttacaagtat
tgaatggttcatttaatatgaaggctgtgtaatgactacagtgaagaagaccaaaaacaa
attcaatattctttggagtctagtaccaccagaagtgactgcttccctccatcaccttat
ttccatgtatttctgaagaaaacctttacttccaacaaagtatttgtcaaatcaaataaa
cctaatcttctgatttaaatttaccatcgagatccaagcttggaatggaaataagagact
atatttgcaaggtaaaacatctcaatcaatttgttaagtagatttggccctagctaaaat
tgagacaaaataaaatcatgccattgttgaaactgttgtg

If you BLAT it against the mouse genome using the UCSC BLAT tool, you will see many hits, most of them almost perfect matches and some of them overlapping the regions given in your bed file:

Your bed file matches the following hits in the image above:

NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3441731   3443690   1960
NM_001085545.1  1960     1  1960  1960   100.0%  chrX   +     3700233   3702192   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3750926   3752885   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   +    31336369  31338328   1960

I don't know why your bed file has all of these various matches for NM_001085545.1, since you haven't told us what the bed file is supposed to be doing. A bed file is just a collection of regions. For all I know, that bed file is just "a few hits for sequence NM_001085545.1 in the genome". There is no reason to assume that each region in a bed file corresponds to a gene.

So there's nothing wrong with your file, you're just dealing with a non-unique sequence from the mouse genome. And you don't have multiple sequences with the same accession, you have one sequence, one accession, but multiple hits for that sequence in the genome.


Is there a standard k-mer count file format?

Not as far as I am aware. The Ray assembler used to (and possibly still does) store the kmers as FASTA files where the header was the count of the sequence, which I thought was a pretty neat bastardisation of the FASTA file format. It looks like this format is also used by Jellyfish when reporting kmer frequencies by the dump command (but its default output format is a custom binary format):

  The dump subcommand outputs a list of all the k-mers in the file associated with their count. By default, the output is in FASTA format, where the header line contains the count of the k-mer and the sequence part is the sequence of the k-mer. This format has the advantage that the output contains the sequence of k-mers and can be directly fed into another program expecting the very common FASTA format. A more convenient column format (for human beings) is selected with the -c  switch.

Jellyfish changed their internal format between v1 and v2 (both not FASTA), because they changed to doing counts based on bloom filters. Jellyfish2 has an optional two-pass method that sets up a bloom filter intermediate file to record kmers, and multiple different final reporting formats.

Khmer also uses bloom filters, but in a slightly different way. It also has been extended to be useful for partitioning and comparing datasets.


Single-cell sequencing dataset has too many barcodes

The BAM file you downloaded has all of the detected 10X barcodes instead of those that represent cells.  In a given 10X experiment the number of input barcodes vastly outnumbers the input cell count.  After the reaction some of barcoded beads are invariably sequenced because they are encapsulated with ambient RNA, dead cells etc. After processing the BAM file Cell Ranger infers which cells likely represent cells and generates gene expression matrices from these files.

You would need to download one of the output files that contains the list of passed barcodes such as the Gene / cell matrix (filtered) and then use the list of passed barocode sequences to filter your BAM file. Alternatively you could approximately select the same cells by selecting the top 2,000 barcodes with the highest numbers of reads.


For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are many papers on this topic.

Error correction alone won't fix all sequencing errors. Remaining errors may lead to bubbles and tips in the overlap/de Bruijn graph. You can identify them by inspecting the local topology of assembly graph.


How to understand and analyse RNA-seq data (for a beginner)?

SRA files are raw sequencing data that have been converted to the NCBI storage format. These can be converted to fastq files, e.g. with fastq-dump from the SRAtoolkit.

Beyond that I suggest you start with the Hitchhiker's GUide to RNA-seq and the DESeq2 workflow to get a background. Be sure to extensively browse the web for vignettes, guides and tutorials. YOu will find helpful resources to get started.


How to highlight specific cells in Seurat 2.4

You can "highlight" any cells you'd like by creating a new column (or modifying an existing column) in the meta data slot within the Seurat object. From there you can modify the contents of that vector to your liking and plot the results.

It's not clear if you've been successful in obtaining the barcodes for the subpopulation to begin with, because that's a separate issue all together, and you need to clarify what attribute you are trying to subset the cells on because there are a number of ways to accomplish that. If you've already obtained the barcodes then I'll describe how you would plot the results below.

Let's assume you'd like to keep the original cluster IDs for all other cells. We'll start off by duplicating the cluster ID column (this isn't necessary, but it's good practice to preserve the original IDs).

n <- seurat_obj@meta.data # for convenience
n$new_column <- n$res.1.0

Now you need to modify the ids for cluster 5 to break them up into 2 populations.

n[rownames(n) %in% my_barcodes, "new_column"] <- "5b"

You can plot the new column using DimPlot with the argument group.by = "new_column", which will make the subpopulation stand out from the rest of the cells in cluster 5.


Retrieving a list of human genes having GO associations

Here's an example for the mouse genome:

library(org.Mm.eg.db)
select(org.Mm.eg.db, c("GO:0048406"), c("GENENAME","SYMBOL"), c("GO"))

You get output like:

                                                    GENENAME SYMBOL
1                                     pregnancy zone protein    Pzp
2 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
3 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
4                                             neurotrophin 3   Ntf3
5             neurotrophic tyrosine kinase, receptor, type 1  Ntrk1
6            furin (paired basic amino acid cleaving enzyme)  Furin
7              proprotein convertase subtilisin/kexin type 6  Pcsk6
8                                                 sortilin 1  Sort1
9                                      alpha-2-macroglobulin    A2m

Edit: Given your comment that you want instead a list of all genes with a GO term associated with them (for human):

library(org.Hs.eg.db) # Install it from bioconductor yourself
unlist(as.list(org.Hs.egSYMBOL)[mappedkeys(org.Hs.egGO)])

The mappedkey() part gets the index into the bimap of all entries with a valid GO mapping. That's then used to subset the symbol bimap (after converting to a list).


They created 10kb bins, which could be called "dividing the genome". These are just 10kb long contiguous stretches (so position 1-10000,then 10001-20000, etc.). This is, of course, a nonsensical way of going about things, because "counting those up" in the regions you posted is the same as taking the width of the regions, dividing by 10000, and rounding a bit.

Presumably you have your markers in something like BED format, so just use bedtools to intersect them and do the counting for you.


Is there a convention or standard of the highest bits in a sequence logo?

This is explained in the original paper [1] describing the logos (emphasis mine):

  The importance of a particular position in a binding site is more clearly
  and consistently given by the information required to describe the pattern
  there [7, 8]. To choose one symbol or state from two equally likely possibilities
  requires one "bit" of information. For example, to communicate the result of
  a coin-ip to someone requires 1 bit of information because only one yes-no
  question needs to be answered: "Is it heads?". If a position in a binding
  site always contains one base (e.g. G), then we need exactly two bits of
  information since two yes-no questions need to be answered: "Is it A or G?"
  (i.e. is it a purine?) and "Is it A or C?". (If the answers to both questions
  are "no", it must be a T.) Furthermore, if a position contains two bases (e.g.
  sometimes A and sometimes G), only one question suffices since a two out of
  four choice is equivalent to a one out of two choice. Therefore, only one bit
  is needed to describe a position in a binding site that contains only purines,
  but two bits are needed to describe a position that always contains adenine.

So yes, 2bits is the highest value you will see on the y-axis of a DNA sequence logo (proteins are a different story since there are more possibilities). In the one you show, the scale is from 0 to 0.2, but it's the same principle.

However, you should think of logos as a useful visual aid, and not a way of getting rigorous mathematical information about a sequence. So don't put too much weight on the bit values. If you need that sort of data, calculate it from the alignment directly.

[1] Schneider TD, Stephens RM., Sequence logos: a new way to display consensus sequences., Nucleic Acids Res. 1990 Oct 25;18(20):6097-100.


Changing a wide range of colours to a limited gradient

TSNEPlot()

TSNEPlot() will always treat your variables as discrete. My approach is to manually generate a gradient with unique colors for each factor level and pass it to the cols.use argument in TSNEPlot().

#generate values for testing purposes, one value for each cell
value <- sample(seq(from=8, to=48, by=1), size = length(rownames(unfiltered_cca@meta.data)), replace=TRUE)
names(value) <- rownames(exp@meta.data)
exp <- AddMetaData(object=exp, metadata=value, col.name="value")
exp <-  SetIdent(exp, ident.use=exp@meta.data$value)
#create a gradient between specified colors, multiply by sequence to get appropriate length
#length.out is the number of different colors, number of factor levels
cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))

TSNEPlot(exp, colors.use=cols)

I would present this with no legend:

TSNEPlot(exp, colors.use=cols, no.legend = TRUE)

FeaturePlot()

You can also simply use FeaturePlot() instead of TSNEPlot() to visualize the gradient. Using the same data as above:

FeaturePlot(object = exp, features.plot = "value", reduction.use = "tsne", no.legend = FALSE, cols.use = c("beige", "red"))

You ask for a continuous scale, but this is not what is shown in your second plot.  You have a different color for each discrete range of values. I assume here that you want a continuous gradient. To achieve something like the second plot you provided requires a different approach.


What are the applications of DNA or RNA pattern matching?

One useful application of pattern matching is for DNA sequence motifs. As an example of a direct application, there is a class of proteins called transcription factors (TFs), which bind to DNA and regulate the expression of genes. There are many TFs. Each TF likes to bind to its particular conserved sequence of DNA, with some wiggle room for some parts of the sequence.
Algorithms for finding and scoring these patterns (and mismatches in patterns) is important for identifying regions where regulation occurs, which proteins are involved, and for answering biological questions related to gene regulation and basic biological functionality — and malfunction, in case of disease.


Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python)

I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:

allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000),
                   ncol=17, byrow=TRUE)
allele1<-allele_counts[ c(TRUE,FALSE), ]
allele2<-allele_counts[ c(FALSE,TRUE), ]
total_counts <- allele1+allele2
allele_counts[c(TRUE,FALSE), ] <-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] <-allele2/total_counts

Edit

I decided to a head-to-head test of this solution vs. bli's pandas solution and your awk solution

I used readr for file IO:

library(readr)
allele_counts<-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)

Here are the timings:

cat Benchmarks/R.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28

The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.


Motif Discovery in Bipartite networks using Cytoscape

I have recently develop an R package (sorry not a Cytoscape solution) to deal with this kind of situations of something grouping other elements, like pathways and genes.
The package is based on GSEABase of Bioconductor to store this relationships. Basically you should store a gene set for each TF and genes it is linked. I would recommend to use the name of the TF as the name of the gene set and then create a GeneSetCollection to be able to use the package.

The package provides with functionalities that allow you to calculate the distribution of genes per gene set (using genesPerPathway), as well as, pathwaysPerGene to know in how many other TF regulate each gene.

I am not sure of the biological inference you can take from this (I designed the package to test this) but you can also see if there is some specific preference for the TF to link to genes which are also regulated by other TF. This could be done by looking at the conditional probability (condPerPathways or sizesPerPathway).
Sorry the manual pages don't have examples and the vignettes are incomplete yet.

Also you could use the package tidygraph to calculate the centrality of the TF or of the genes you want.


Possible?  Sure.  Your biggest problem is that 10x is expecting barcodes from its whitelist.  I used to alter the python scripts to accept a 'new' kind of chemistry where the barcodes and umi lengths matched what I had, but I found that my modifications stopped working when the version updated,  and I couldn't figure it out, so I started modifying fastqs instead.


Given a .gb file and a locus - how to get relevant annotations in Python?

Biopython is the main package for this. It's only a few lines, so it is not reinventing the wheel. Unfortunately, this feels like a homework question, so nobody can write the code for you as per no homework policy.

But pointers are okay...

So you want to iterate across the features table (a list in Biopython) of a record and find the case where the qualifier['locus'][0] matches your query.

The things to watch out for are:

filter by type (CDS?)
entries where there is not key among the qualifiers, so add error catching (try: ... except Exception: pass
the values of the qualifiers of a feature are a list. so add a [0]
multithreading.Pool might help if speed is a worry (which it really shouldn't). Asyncio and thread are not the way (both single core).


Is Wikidata a reliable KB for genomic and drugs annotations?

The question is what you mean with "reliable". For bacterial genome annotations there is one Wikidata based platform called WikiGenomes (publications).


Splitting a GenBank file into smaller files

You can use NCBI EntrezDirect to download data in GenBank format for a specific region of the sequence as follows:
efetch -db nuccore -id NC_000962.3 -format gb -seq_start 1 -seq_stop 99000 > seq.gb


Creating frequency plot from weblogo with frequency unit on the Y axis

Within weblogo, add the --units probability argument and edit the figure with Illustrator (PDF/SVG) or Inkscape (SVG) to modify the y-axis, changing fractional probability ([0,1]) to a percentage ([0,100]) by multiplying by 100. Refer to the manual for more detail.


How to verify if the stationary states of a biological model are correct?

After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.

But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct.

Regarding the model is taking forever, I took the logistics part and ran well!


What is the significance of enrichment factor regarding machine learning methods?

I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is...

Where possible, it's best to use explanations provided in the original source. This is because definitions could be different, depending on who is writing. Here's some descriptive information from the paper that may be more informative:

  The enrichment factor (EF) is another metric we use and is very
  popular for use in evaluation of virtual screening performance. EF
  denotes the quotient of true actives among a subset of predicted
  actives and the overall fraction of actives.... To evaluate the
  maximum achievable enrichment factor (EFmax), the maximum number of
  actives among the percentage of tested compounds is divided by the
  fraction of actives in the entire dataset. To more consistently
  quantify enrichment at 0.1% of all compounds tested, we report the
  ratio of EF/ EFmax. This normalized metric is useful because EF values
  are not directly comparable across different datasets due to the
  maximum possible EF being constrained by the ratio of total to active
  compounds and the evaluated fraction as shown in the equation above.

So for the purpose of this paper, they're saying that this ratio is an appropriately normalised substitute for a plain enrichment factor. High enrichment factor == good; low enrichment factor == bad.

Given the information you've provided, I'd agree that the curves in this figure appear to suggest that the multitask DNN performs worse with smaller numbers of active compounds (where Naive Bayes or Logistic regression would be a better choice). It's hard to tell from the graph (and legend) what the curves represent, and the clipping at EF/EFmax = 0 doesn't help.

The curves have a strange shape, suggesting the number of active compounds is dependent on the EF/EFmax ratio, rather than the other way round (which would make more sense to me). If I consider the curves to be something like a smoothed mean of the active compound number, then any place where there are fewer dots (i.e. low EF/EFmax for Multitask DNN) will have more associated error. For Multitask DNN, I see a large gap between 0.0 and 0.2, with lots of outliers on the 0.0 line (for all methods), so presumably the 0.0 outliers are skewing the curves.

Given that this is a preprint, the generation of the curves isn't specified, and the other figures (especially Figure 3) seem to tell a different story (i.e. that multi-task DNN performs better at all scales), I'd recommend contacting the authors to ask / clarify what the curves represent, or ignoring the curves entirely.


Accidental mapping of eukaryotic reads in a metagenomic dataset

If you know the eukaryotic contaminant present you could use bbsplit.sh from BBMap suite to split/bin reads first using a reference for that contaminant (into one or as many bins as reference sequences you provide).


Model the effect of treatment on mutation burden

Depends how you consider "linear regression", if you mean Pearsons, Spearman's analysis then general linear modelling (GLM) is a far better approach if you have continuous data. I don't really understand TMB, but that is a critical output then GLM will work. Logistic regression is an alternative approach data type depending. You have to perform a Q-Q plot and the residual against against any transformations you make (you might need to transform TMB). The lower the residual the better, the straight the Q-Q plot the better.

To explain how GLM works long hair and height correlate in a straight regression analysis. The implication is if you want to be tall - cut your hair! Of once you account for gender and send it through a GLM then gender and hairlength will correlate far stronger than hair length and height.

Mac OSX has a very easy to use pointy clicky, but highly professional, GLM package on its AppStore called Wizzard (bit expensive). It doesn't do general linear mixed modelling however. Everything does GLM of course, its just this package is really easy to use and you can figure out GLM by trial and error.


Assuming none of the genes have duplicate entries across chromosomes (or strand or other biologically implausible things) in R one can:

library(rtracklayer)
library(GenomicRanges)
gtf = import.gff("file.gtf")
grl = split(gtf, gtf$gene_id)  # This produces a GRangesList
grl = endoapply(grl, function(x) {
    foo = x[1]
    foo$type = 'gene'
    start(foo) = min(start(x))
    end(foo) = max(end(x))
    return(c(foo, x))
    })
gr = unlist(gr)

You can then either write gr out to a file (either with export() or just printing the rows).


Interpreting MultiQC plot of GATK BaseRecalibration data?

Turns out two recalibration tables are needed to show pre-recalibration and post-recalibration plots. I previously only had generated the pre-calibration table by running BaseRecalibrator on the BAM that hadn't undergone ApplyBQSR to adjust base scores. Here are the steps to show pre and post recalibration plots in MultiQC:

Run BaseRecalibrator on your initial input BAM. This will create a recalibration table with the original base scores.
Run ApplyBQSR with the recalibration table from step 1 to create a new BAM with the adjusted base quality scores.
Run BaseRecalibrator on the adjusted BAM generated in step 2 to create a new recalibration table with the adjusted base quality scores.
Use both recalibration tables as inputs to MultiQC. This should create both plots.


HTSeq-count is returning 0 for every gene, instead of expression value

Maybe the contig names in your BAM files don't match those in the GFF file? Ensembl uses contig names like 1, some other reference sources use chr1.

Generally, for ensembl GFF files, HTSeq's default value for idattr should be fine though.


Virus genome: protein inside another protein

Note that the protein sequence of ORF8b is not contained within the protein sequence of protein N. That means that they are translated from different frames (though they are in the same orientation). Indeed, if we just subtract the start coordinates we do not find a multiple of 3.
Thus, they cannot be results of the same translational event. This means that the smaller protein could result from translational slippage (as suggested by Matteo Ferla), "leaky scanning" of the translation machines for start codons, or it could simply result from a different ribosome assembly site, recruited via the viral mechanism of an internal ribosome entry site, which allows a ribosome to be recruited in a mRNA cap-independent fashion (though this may be less likely).
For some mechanistic details of how one such ORF is translated, see this paper.
It looks like such things are sometimes called "subgenomic mRNAs". For a paper on such nested genes in coronavirus, including I believe Orf8b, see here.


Making PLINK compatible files from VCF file without phenotype information

You can use --pheno to specify a phenotype file you want to use, the 1st and 2nd columns being FID and IID. These are just used on the fly and cannot be written to the BED files as I understand. --make-pheno is used to write binary phenotypes, which is unsuitable if your phenotype is actually continuous.
But I still do not quite understand what you want to do here.


Alignment to verify plasmid from forward and reverse Sanger sequencing reads

Align all reads from a sample at the same time IMO. This will give you all of the information at once, and it should not be difficult. It should be very obvious whenever there's a disagreement between reads or between reads and reference.
I would also recommend inspecting chromatograms directly at all sites of interest (such disagreements, insertion sites, possible errors or mutations).


How low does an E-Value need to be for a gene to code for the same protein as query?

I assume that you mean "is homologous to" by "encodes the same protein as". Identical ("same") proteins should be fairly trivial to identify using a tool such as tblastx.
You should probably not use E-value for this. E-value is similar to a p-value, in that it represents some aspect of the data (in this case a sequence and a database) compared to some null hypothesis (meaning in this case statistics on the database). Therefore, it depends both on the specifics of the data and the hypothesis, and there is no universal threshold that will work across all sequences and all databases. (See discussion of the meaning of e-value.)
E-value does not in any way represent the content of the sequence or its relationship to a sequence in the database. You could have E-values close to zero or E-values > 10 for a perfect sequence match or a very loose sequence match, depending on the input sequence and the database.
However, BLAST output does include other information that might be helpful. You could (for example) look at the length of the alignment and the number of mismatches/identities in the alignment. This will give you a pretty good guess at how similar the two sequences are; unless they are very divergent you would expect them to align across most of the length with somewhat high protein identity or similarity (>25% probably?). This process will still not work 100% of the time, but I think it is closer to your intent than any given approach involving E-values.


What is the recognized way to reference accession data (.sra files) from NCBI as a URI?

The nomenclature is no longer being followed as NCBI is moving all the files to S3/Google cloud. pysradb allows you to fetch the metadata and get the URLs for each individual SRR. See examples in this notebook.

You can also download an entire project by doing:

pysradb download -p <SRP>


scRNA-seq multi-dataset integration for small datasets

TLDR: Harmony can work on 106~ cell samples but has less versatility then methods like BATMAN. BATMAN is useful if you need your data for differential gene expression and single-cell eQTL and can work on ~200 cell data.
Well their certainly is a minimum which comes down to several factors. As many of these rely on algorithms like KNN and other clustering algorithms for integration the threshold is dependent on the quality of the data and the  algorithms abilities to predict reasonable clusters from that data. As well as the downstream uses of the data...
BATMAN highest performance on simulated set 200 cells and 1000 cells
With 200 cells from one batch and and 1000 in the other only BATMAN succeeded in clustering the high dimensionality simulated data well. If you need your data for differential gene expression and single-cell eQTL, BATMAN seems to be the correct choice, for small sets. "The downside of these methods is that they operate in latent space, which limits their interpretability and use in downstream analyses such as differential gene expression and single-cell eQTL analyses."
"A shows that the two original datasets become better mixed with each other. BATMAN has the highest performance not only when considering the top two principal components but also when considering more top principal components. It is the only method that manages to efficiently maintain a high iLISI score in a larger number of dimensions."
Harmony is best if... differential gene expression and single-cell eQTL is not necessary
"Moreover, we show that Harmony requires dramatically fewer computational resources. It is the only available algorithm that makes the integration of ~106 cells feasible on a personal computer."
Computing time benchmarks-section
It seems as the runtime of these algorithms is less than 2 minutes, using a set with less than ~2000 cells will not benefit your computational time much.
"When processing a small data set of ~ 2000 cells, all four methods took less than 2 min."
"To obtain such datasets, we downsampled the MCA and TM datasets to obtain a total of 9 sets of data containing between ~ 2000 and ~ 140,000 cells, while the number of highly variable genes (HVGs) was controlled in a range from ~ 2000 to ~ 3000 (Table S1)."


Using cellranger mkgtf w/ NCBI GFF3 files

Regarding your comment:

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

The link you have provided corresponds to an entry in Genbank format and not to GFF format. I believe you need a file in GTF format in particular.

Nevertheless, given the fact that you have manually added the corresponding entries to the reference fasta and gtf files prior to mkref, you should have seen your gene of interest in the output. I would double check one of the original fasta headers and original gtf entries in the reference file and match these for the new entry.


choosing a good representative genome subset

Check out the progenomes database: http://progenomes.embl.de/


How to map gene names to ensembl Ids in R

You can pull a full list of Ensembl gene IDs and MGI symbols from biomaRt in R. From there it should be easy enough to filter for the genes you want:

library(biomaRt)
hsapiens_genes <- getBM(attributes = c("ensembl_gene_id",
                                       "hgnc_symbol"),
                        mart = useMart("ensembl", dataset = "hsapiens_gene_ensembl"))

You do not find some genes in your example because you were not using the official gene name. For example CTGF is offocially called CCN2 and with this name you will find the appropriate Ensembl gene ID in the biomaRt table. This is why I personally prefer to use gene name that are like ENSG00000118523|CCN2 in my work. That way you always have a unique identifier (the Ensembl ID) that is unambiguous and you do not run into trouble that for one gene many more or less official gene names exist.


How can I extract all known mutations of my BAM (or SNP/INDL files)?

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


Difference between copy number neutral reads and active reads

In genetics normally one gets a copy (allele) of a gene (or actually a chromosome) from each parent. If the alleles differ between parents, for example by a SNP, one would have one copy of each allele. With (NGS) sequencing you would get an even amount of reads for both copies. However with somatic mutations, hence mutations which occur after meiosis and gamete forming, one can lose or gain parts of the chromosome (copy number variation). In this situation, some genomic regions might be deleted or amplified. So the allele of one parent can be deleted, leaving you with only the allele of the other parent.

As requested, I have made a simple drawing of a situation where a region differs between parent alleles with SNPs. After a deletion of one of the allele only reads of the other allele are left.


How to show bootstrap values on a phylogenetic tree constructed with RAxML

Depending on the search parameter settings RAxML produces several output files among which you normally find:

  RAxML_bipartitions.result: If you used the -f b option, this file
  will contain the input tree with confidence values from 0 to 100 drawn
  on its nodes! It is also printed when -f a x have been specified, at
  the end of the analysis the program will draw the BS support values on
  the best tree found during the ML search.

  RAxML_bipartitionsBranchLabels.result: Contains the same information
  as the file above, but support values are correctly displayed as
  Newick branch labels and not node labels! Support values always refer
  to branches/splits of trees and never to nodes of the tree. Note that,
  some tree viewers have problems displaying branch labels, they are
  however part of the standard Newick format.

from https://cme.h-its.org/exelixis/resource/download/NewManual.pdf

If you don't find the files in your output directory you did not specify the above parameters. If a rerun is not feasible you could you use the multinewick file RAxML_bootstrap.result with raxml or as input for the program splitstree to generate a consenus tree with different rules (majority rule, strict, etc.). However, these options don't guarantee to generate the topology of the maximum likelihood tree.


Pointing Cromwell to local EBS mounted storage on AWS Batch

Somebody helped me with this one through other channels, but I wanted to make it publicly available:
Make a tarball of the STAR index. The option --strip-components=1 strips the first level of the path from the untarred items; thus, the expected tarball was made with something like
tar -cvf ${star_index} the_star_index_directory/

Then point Cromwell/WDL to the tarball, e.g. on S3. The GTEx container will unwrap the tarball on the fly during execution.


Tools to create annotated table of variants from VCF

Hail might be an option for you.

It is actively developed by a growing team at the Broad.  It is rigorously tested (continuous integration, continuous deployment, bug reports get regression tests, blah blah blah).

It was designed to solve this problem (among others).  It can import a variety of formats, including VCF, TSV, UCSC BED, JSON and interval files.  (We don't have explicit support for GFF, but we can probably handle them with general facilities.  If not, get in touch and we'll add support.)  It can call out to VEP (and soon Nirvana, Illumina's VEP rewrite).  It has general facilities to transform, filter, clean and query data.

What's more, we've curated a large collection of annotation resources (currently ~22 databases of annotations + VEP), hosted in a public bucket in Google cloud, and built an interactive query builder to select which resources you want to use.  Get in touch if you'd like us to add additional resources.  You can read more about it here.

Disclaimer: I work on Hail.


Given a VCF of a human genome, how do I assess the quality against known SNVs?

To achieve (at least some of) your goals, I would recommend the Variant Effect Predictor (VEP). It is a flexible tool that provides several types of annotations on an input .vcf file.  I agree that ExAC is the de facto gold standard catalog for human genetic variation in coding regions.  To see the frequency distribution of variants by global subpopulation make sure "ExAC allele frequencies" is checked in addition to the 1000 genomes.

Output in the web-browser:

If you download the annotated .vcf, frequencies will be in the INFO field:

##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|SYMBOL_SOURCE|HGNC_ID|TSL|SIFT|PolyPhen|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|AA_AF|EA_AF|ExAC_AF|ExAC_Adj_AF|ExAC_AFR_AF|ExAC_AMR_AF|ExAC_EAS_AF|ExAC_FIN_AF|ExAC_NFE_AF|ExAC_OTH_AF|ExAC_SAS_AF|CLIN_SIG|SOMATIC|PHENO|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE

The previously mentioned Annovar can also annotate with ExAC allele frequencies.  Finally, should mention the newest whole-genome resource, gnomAD.


Are there any resources for calculating / representing evolutionary relationships within gene clusters?

Firstly you wanted a webserver to build a tree, ok, I would use PhyML server. The good thing is that with only 22 taxa the calculations are easily doable.

What I'm assuming is you can map your rearrangements to produce a single homologous gene alignment for all taxa. So you appear to infer you can break your 20-30kb contig into individual aligments of homologous genes which has incorporated the rearrangement and enables easy production of an alignment file. This would leave you with around 20-30 gene alignment files. I hope that is correct, because that is how you get around the rearrangement problem and enables the data to be reassembled for the final tree (below).

Doing the tree is simple, because you can exploit the assumptions used in tree building regarding independence between mutations at each alignment position. The theoretical proof, the data is sound, is complicated. In summary,

build a consensus phylogeny (one tree representative of lots of trees) using two different approaches
combine the data into a single alignment, if the above looks good, or excluding the alignments that are affected recombination.

The theoretical issues with what you want to do are complicated because you need to distinguish bifurcating (clonal) homologues from recombination that could have occured in the reassortment events. The following is a  purists approach, you might not want to do this stuff, but it gives you some idea.

Just before I do that ... I just want to check the alignments assumption is correct? and which bacteria?


What methods are available to find a cutoff value for non-expressed genes in RNA-seq?

I'd like to find genes that were not expressed in a group of samples and were expressed in another group.

This is, fundamentally, a differential expression analysis, with a twist. To solve this, you’d first use a differential expression library of your choice (e.g. DESeq2) and perform a one-tailed test of differential expression.

Briefly, you’d perform the normal setup and then use

results(dds, altHypothesis = 'greater')

To perform a one-tailed test. This will give you only those genes that are significantly upregulated in one group. Check chapter 3.9 of the vignette for details.

Of course this won’t tell you that the genes are unexpressed in the other group. Unfortunately I don’t know of a good value to threshold the results; I would start by plotting a histogram of the (variance stabilised) expression values in your first group, and then visually choose an expression threshold that cleanly separates genes that are clearly expressed from zeros:

vst_counts = assay(vst(dds))
dens = density(vst_counts[, replicate])
plot(dens, log = 'y')

(This merges the replicates in the group, which should be fine.)

Counts follow a multimodal distribution, with one mode for unexpressed and one or more for expressed genes. The expression threshold can be set somewhere between the clearly unexpressed and expressed peaks:

Here I used identify(dens) to identify the threshold interactively but you could also use an analytical method:

threshold = identify(dens)
quantile = sum(dens$x < dens$x[threshold]) / length(dens$x)

# Using just one replicate here; more robust would be to use a mean value.
nonzero_counts = counts(dds, normalized = TRUE)[, replicates[1]]
nonzero_counts = nonzero_counts[nonzero_counts > 0]

(expression_threshold = quantile(nonzero_counts, probs = quantile))

26.5625%
4.112033


How do I estimate branch length of a maximum parsimony phylogeny in R?

@Nathan S. Watson-Haigh's answer does give branch lengths, but NOT parsimony branch lengths, which may be what the OP wanted.  nnls.phylo is giving distances calculated from a pairwise distance matrix. To get the parsimony branch lengths (i.e. number of changes mapped on to each branch), score the tree against the dataset using acctran:

treeRatchet <- pratchet(primates)
treeRatchet.BL <- acctran(treeRatchet,primates)
plot(treeRatchet.BL,main='Phylogeny (Parsimony branch lengths)')

Note that the branch lengths in the parsimony phylogram are 'chunkier' than those in the nnls.phylo phylogram -- that's because

they are all integers, and
they aren't corrected with reference to a model of substitution that considers multiple hits or base composition.


Transform traditional blast output to `--outfmt 6`

The blast_formatter command allows you to produce output in any of the formats supported by BLAST. However, it requires that you run blastx with --outfmt 11 (ASN.1 format). I'm not sure there are any tools that will do the conversion from the full default BLAST report. This could potentially be scripted, but would be complex enough that you're probably better off re-running BLAST.


Differences in webtool and standalone Clustal Omega Percent Identity Matrix

Try clustal omega, instead of clustalw because it is really old. It should report to 2 dp. Variations do occur depending on how indels are treated, but this is just rounding.

You could try

alignment -tree -pim -score

The idea is to convert to absolute rather than percentage and might overcome the rounding error


Calculating the mortality rate of a pandemic, e.g. coronavirus 2019-nCov?

This may not strike most as a bioinformatics, but getting the key clinical outcome is essential in understanding the molecular basis of pathogenicity.

I think the mortality rate is over-reported. This is not to say the situation of 2019-nCov is not serious - it is very serious.

The two essential factors missing in your equation for 2019-nCov are:

Age. Mortality is likely skewed towards the old and very young
Seroprevalence. To identify the number of asymptomatic patients.

Old people have a weakened immune system (e.g. reduced T-cell response). The very young have an underdeveloped immune system. The problem is hospital cases are skewed towards those with severe infection, which is likely the old and very young. In 2019-nCov the first infecteds were middle aged men, but likely to be restricted to those exposed to the "Wuhan meat market".

Problem
Asymptomatic "spreaders" are known in influenza and are suspected in 2019-nCov. Asymptomatics will inflate the mortality rate (because we're not detecting them).

Solution(s)

One solution is to conduct seroprevalence surveys. Serological surveys screening for antibodies (IgG) against 2019-nCov will identify asymptomatic patients, using the first formula you stated. However, would you seriously go around the general population in a 2019-nCov epidemic taking blood samples? Errhhh no.
The bio-statistics approach is to use a catalytic model and critical to this is the age of the patient. Catalytic models are a priori models that incorporates the age-prevelance profile of infection. By introducing the incidence and age of hospitalisation of e.g. the old and the young, the rest of the distribution, i.e. asymptomatics can be calculated by regression. This estimates the total population of infecteds and again you use equation 1 to calculate the mortality rate recorded by the hospital.
Mathematical models based on transmission dynamics. The problem here is the rate of transmission is unknown.

I would assume the hospital is the place virus-induced mortality is reasonably accurate.

The second formula you described works very well for Ebolavirus, no-one has ever heard of asymptomatic Ebolavirus. It would also work well for another betacoronavirus Middle Eastern respiratory syndrome. The key problem with 2019-nCov is we don't know the number of asymptomatic patients or patients with such mild infections they never report illness.

Summary
In summary, I think the mortality is being over-reported if it is based on your equation 2. Equation 1, however, will also over-report the mortality rate.


Most likely, the GWASs that generated your summary statistics used other imputation panels than 1000G, like HRC. Clearly, PLINK can't estimate the LD for SNPs that aren't found in the reference, and most likely just ignores them. I think the easiest solution is to use whatever individual-level dataset you have to estimate the correlations - even if it's small or from a different population. The estimate will be biased, but you don't really need high precision for clumping.

Then again, you can skip PLINK and try searching for SNP prioritization methods based on summary statistics - it's an area of active development, and something might be available to suit your needs.


How can I only get the species name for fasta sequences from blast results?

There might be a "cleaner" way to do this, but what I would do is first remove all of the spaces from the fasta sequence names and replace them with underscores using something like sed.

sed 's/ /_/g'

Run your BLAST after you remove the spaces and you will get the whole ID string in the result instead of only that first part (e.g. NP_567483.1). From here you can parse out the species name, or whatever other ID you would like.

Hope this is helpful! :)


Getting System.IO.FileNotFoundException with MaxQuant running on Linux with mono

...seems the mqpar file was not compatible with the maxquant version.
Make sure the mqpar.xml file shows the same version number as the MaxQuant binary.
MaxQuant is not really made for the command line. So, the errors shown are sometimes quite misleading.


Is it possible to calculate p-values for fold changes of single replicate RNA-seq samples?

DESeq2 is able to adjust p-values for single-replicate samples by estimating shared dispersion across all samples. DESeq2 will give a warning, but try its best to carry out the analysis. More information on doing that can be found in the DESeq2 vignette and the DESeq2 manual page.

The results will not be anywhere near the sensitivity that can be achieved when replicate information is available, but it might be good enough if only very obvious large-scale changes are of importance.


I want to advise you to use FASTA format, instead of obscure unknown formats. In R, there are many packages that can handle FASTA format and it is therefore more convenient. I will show you how to make a FASTA file from your data, and then I will show an example of how to import such a FASTA file and sub select your proteins of interest.

input <- scan("seq2.txt", what = character(), sep = "\t")
length(input)
new <- gsub("f", ">f", input[-1:-2])
write.table(new, "seq2.fasta", quote = F, row.names = F, col.names = F)

library(seqinr)
fasta.file <- read.fasta("seq2.fasta", seqtype = "AA", as.string = TRUE)
prot.id <- scan("protein_id.txt", what = character(), sep = "\t")
selection <- fasta.file[names(fasta.file) %in% prot.id]
length(selection)
[1] 87

head(selection, 2)
$`fig|6666666.151743.peg.5044`
[1] "MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAITSEQAIKLGITQPDVIAQMVPGLQLNHSAWGLTPFLRGVGTNSVAIGNEPAVALFVDDVYLPTGSAGVFQLNSISSIEVLKGPQGTLFGRNATGGVIHIRTRDPSFDPSASATASYGNYDTISGQLYGSTGLTDTVAINFALQGIDQRKGWGVNGTTGNEANSMSGWGGRAKLLWNMGEETSLLVAGSYDRRSSDQGQPGRSIAGTFTRDGYSADAAGLGFYDVMSNLDTVQKSTTRSLSGKLTHDFGSARLVLISAFSKTVVDPVVSDSDGSPSPFFQTREFAMGGRTFTQEVQLLSPEGSAVKWILGAFYLHDRSFMDATYEGQAFGTSKSHTDAVQRTNSIAGFGQVTTEILPRTNLTLGLRYTSDRRSLDAVARLLPATGAAVTFGPFEDAKTFNSVTGRASIDHHFSDDLMGYIAYNRGFKSGIYNFTGIVAGTTTLPPPVNPETLDAYTIGFKSEFLNRRLRVNAEAFLYEYSDLQVQIVLPTGTTLTNASKATIKGIDFEISAIPVRNLTVSASIGLLDGRYKDFRNAPSFFLLPPNAPTPYDPARCSGPNPYPAVQGSAPATQRSCDVSGNRLIQAPPFSSNLQLLYTLPTDLGDFDLSASWAHGGNYEYDVFSSPLTRQPRTDIVNGSVQWTAPGGTFNVRLWANNLTGEKYYTIATTGSTAGAKYAAAAPRTYGVTLGANF"
attr(,"name")
[1] "fig|6666666.151743.peg.5044"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5044"
attr(,"class")
[1] "SeqFastaAA"

$`fig|6666666.151743.peg.5063`
[1] "MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGTKREEPLARLPISVAIVQAPGEQGGALLPGTEEVAQRLDGLTLTNLGPGRNRQFIRGVADSPFNGNSQSTVAILINDARATYNAPDPDLKLVDVERVELLKGPQGSLYGTGTLGGVYHIVTRPPDAGSFAAMAGMGASLATHGEMGGSAEAMVNLPVVADRLALRAVAYAVREPGWIDTAGRGDANRGELRGGRVALRALAGDWTVDLNGMAQMRDVADSQYVYARHALERPAQPAEPHDNDFALMALRVAGKLGGADLVYAGNWVDHDVSSTLDASDLSASAFGIAPPALYLEQRRNRLFDQEVRLSQTRDDGLSWLAGLSFVDVRNHLSGVVRGVRDETVGTGIQGVVEFSAFGDMRMPLTGALSLSVGGRLFHSRIEEEQGARHARIEASVARSGFTPSLALSWQAGAHDFLYARYAGALRPGGLSRTGAASAERFESDELQTVEMGWRHSASAWLTVNLDAFLTTWSHMQSDYLLADGLIGTRNVGDARIRGVEASVDWKPAAGWAISAGLVGQSALLTRTNEGLRLEDRHLPVVPDWKARFGVNRHFDLSGWSGSMGVNLAYLGPSRLSFDPGLDRRMGDYALANADLALRRKDWSIRLSGSNLLDCKADSFAFGNPFALARGTQYTPLRPRQWAIGFTKNW"
attr(,"name")
[1] "fig|6666666.151743.peg.5063"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5063"
attr(,"class")
[1] "SeqFastaAA"


Calculating average coverage for .bam files (sequence data)

For a quick estimate you’re making it more complicated than necessary.

The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size.

samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list.

Putting this together, we get
samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'

(In this case I’ve used 120 as the read length.)
However, note that this over-estimates empirical sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average mapped length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.

If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the -a argument to samtools depth: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:
sum=$(samtools depth "$1" | awk '{sum+=$3} END {print sum}')

Some more remarks on your code:

The -P flag for grep is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: grep '^@SQ' works just as well.
While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:

You can often combine cut and awk, in this case by setting the field separator:
awk -vFS=: '{sum+=$3} END {print sum}'

You can often combine grep and awk:
awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'

Always double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters are valid characters in paths!
Your code is clean and readable. Well done.