scRNA-seq multi-dataset integration for small datasets There have been a few methods proposed for integration (or batch correction) of scRNA-seq datasets, such as Seurat CCA, MNN Correct, Scanorama, and Harmony. The concern is generally about the maximum number of cells that they handle, but I haven't seen any discussion about the minimum number of cells. I am confident they can all handle 10k cells reasonably well and will fail with 10 cells, but where do you draw the line? Is there a method that works best for small datasets?

For example, with plate-based platforms like Fluidigm, many experiments only have 96 cells and potentially much less after quality filtering. How can those be used?

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


State of the art mutation simulation software There are many features affecting mutation probabilities, e.g. CpG mutations are 10-fold more likely than other types of mutations.

Is there a model (preferably with software) which can take two aligned genomic regions, estimate parameters of the neutral mutation procces and then simulate mutations for another region? Obviously, there are models like HKY and GTR, but they are somewhat simplistic and do not take the context into account.

One model I found is this, but it seems that there is no software implementing the model. Is there something more recent?

I am also looking for software which can account for indels.

Question I suspect this question is about the phylogenetics of methylation and the approach the investigator is proposing would be the last approach to use.

Summary
The approaches to assess the phylogenetics of methylation in order of preference are:

dN/dS between CpG sites and non-CpG sites,
An explicit molecular clock between CpG sites and non-CpG sites
A likelihood ration test (LRT) based on a null distribution of mutations generated by a Monte Carlo algorithm

Background/Rationale
The approaches you are using have been around for a long while and generate the random distribution of mutations for a given phylogeny. Its phylogenetics in reverse and the technique is known as Monte Carlo, you start with a randomised probability and send it through a parameterised model to predict the amino acid/nucleotide. Thus it is ML(max likelihood), Bayesian, "phylogenetic-HMM" in reverse. It is used within likelihood ratio test calculations. The model is determined by a standard ML, Bayesian phylogeny algorithm, i.e. its circular mathematics because there is no independent means of calculating the mutation behaviour, so the precise context you are using this for needs to be carefully considered.

Packages There are lots of packages that do this type of Monte Carlo simulation, SOWHat being a good one (LRT) and you generate a large number of replicate data sets (100 or 1000). The one of the base algorithms is "seq-gen", although PAML may have implemented this.

Considerations To use this style of approach you need to carefully consider your question. The "mutation simulations" when analysed via a phylogenetic program will produce the same parameters that you initially set and trees of basically the same length. If you are using this to generate a null distribution for a phylogenetics test these approaches are useful, you then compare the observed likelihood against the null distribution. If you are using it to work out whether the mutation rate at CpG sites is higher than at other sites, its one approach amongst a number of alternatives.

Drawbacks Calculating a de novo null distribution of the mutation rate is computationally very expensive and therefore tends to the last calculation performed. It will address a single hypothesis, usually regarding topology.

Crux of the problem To really do what (I think) you want to do you need an independent measure of mutation and mutational behaviour and its not trivial to achieve. You would have to consider it like a machine learning calculation, with a formal training set, I don't know whether such an approach has been implemented.

Strengths/Summary In summary, Monte Carlo phylogenetics simulation has a limited application because of circularity of the calculation, BUT, BUT, BUT what I have omitted is when it is appropriate it is a very powerful test indeed.

Points of interest I do like the population dynamic simulations within Beast, this certainly "a priori" and its potential IMO needs further exploration. However, I don't think you looking at how molecular epidemiology might impact the mutational behaviour methylation. I have not looked at Indelible, but looks interesting.

I was going to talk about "processor war simulations" but its probably off-topic.


Why DoHeatmap Does not show all genes in genes.use? I am heatmaping a list of genes by DoHeatmap function in Seurat R package. I am sure I have 212 genes but heat map shows only a few of my genes

> DoHeatmap(
+        object = seurat,
+        genes.use = genes,
+        slim.col.label = TRUE,
+        remove.key = TRUE)
> dim(as.matrix(seurat@data))
[1] 12293   209

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

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

Please look at heat map, I was trying to plot 212 genes but you seeing only few genes are among my 212 genes

Please somebody save me from this confusion

https://github.com/satijalab/seurat/issues/287

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


How to calculate lipophilicity of proteins? The question says it all. I want to calculate lipophilicity of proteins to screen the compounds out for antigen-antibody interactions.

Oligopeptide

The lipophilicity of an oligopeptide or a small molecule can be calculated with RDKit, the classic compChem python toolkit. The measure is logP or partition coefficient, namely the magnitude of the ratio of the concentration in octanol over than in water, so the higher it is, the more it partitions in the octanol phase. However, with a compound of several dalton it may be good to see what the contribution of each atom to this score are. This measure is the Crippen contribution to logP, which can be handily plotted on the molecule. Note that these are may be at first glance look similar to Gastaiger partial changes, but are not.

Okay, first the imports

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from rdkit.Chem.Draw import SimilarityMaps

Now let's make the molecules. Let's start with something we can simply use a SMILES string of the web: Bortezomib.

name = 'Bortezomib'
mol = Chem.MolFromSmiles('O=C(N[C@H](C(=O)N[C@H](B(O)O)CC(C)C)Cc1ccccc1)c2nccnc2')

Now we can do all the relevant calculations on the mol object:

print('logP', Descriptors.MolLogP(mol))
AllChem.Compute2DCoords(mol)
contribs = Chem.rdMolDescriptors._CalcCrippenContribs(mol)
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, [x for x,y in contribs], colorMap='BuPu', contourLines=10)
fig.savefig(name+'_crippen.png', bbox_inches='tight')

LogP=0.3

colorMap='BuPu' is the matplotlib colorscheme. I have no idea which work best, I have spent way too much time trying and I do not know. But with BuPu, white is hydrophilic, purple is hydrophobic.

Say the oligopeptide does not have a SMILES or we don't want to do it that way. Let's take β-Casomorphins 1–3:

name = 'β-Casomorphins 1–3'
mol = Chem.MolFromFASTA('YPF')
mol = AllChem.ReplaceSubstructs(mol,
                                Chem.MolFromSmiles('N'),
                                Chem.MolFromSmiles('[NH][OH]'),
                                replacementConnectionPoint=0)[0]

The result with the last block from before gives:

logP=1

Proteins

With protein, it is a lot harder. Membrane affinity of protein is dependent on their structure. In PyMOL you can use the Adaptive Poisson-Boltzmann Solver (APBS) to show the surface by charge. For transmembrane protein there is clear ring of hydrophobicity —CHARMM-GUI membrane builder and a variety of tools for example can guess where it is. For a membrane associated protein only one side will be hydrophobic with a tell-tell charged ring around it (e.g. DHR1 domain) and there are no hard and tested methods for that.


GWAS MAC filter Interpretation I am performing a GWAS analysis and try to understand the influence of the minor allele count filter.Setting the filter to 1 % gave me this plot and I am confused about the same -log10 pvalues around ~ 3.8 which shows a strange pattern (many dots in one line with same value).

If I increase the filter to 6 % the dots in a line disappear but I loose also some significant SNPs:

Sample are about 77000 SNP positions (barley), unimputed from GWAS (all NaNs removed). GWAS was performed with 200 phenotype samples.
How to interpret the dots in  a line?

This looks like genotyping failed for some samples. It's unclear from the op whether the genotyping data was controlled for quality, but I would recommend performing sample QC.
On a side note: in human gwas with imputed snps, I would expect to see clusters of significant snps as a sign that the significant snps are not an artifact. In the op there are some single significant snps. This might be a sign that those significant snps are artifacts. There are also some snps that form clusters and therefore look more promising.


Dmpfold error "No such file or directory: 'best_qdope.pdb'" when running "run_dmpfold.sh"? I am running DMPfold, "seq2maps.sh" ran successfully, but whrn running "run_dmpfold.sh", it gives me this error:

   Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
seed =  30603
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/run_dmpfold.sh: line 177: c600m2-tail-before-lysin-sequence_[1-9]*.pdb: No such file or directory
No atoms read from PDB file
0 models read from PDB file
Outliers:

Traceback (most recent call last):
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-softmax/pytorch_dmp_iterdistpred_cb.py", line 82, in <module>
    main()
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-softmax/pytorch_dmp_iterdistpred_cb.py", line 55, in main
    with open(sys.argv[3], 'r') as refpdbfile:
FileNotFoundError: [Errno 2] No such file or directory: 'best_qdope.pdb'
Traceback (most recent call last):
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-hb/pytorch_iterdmp_hbpred.py", line 75, in <module>
    main()
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-hb/pytorch_iterdmp_hbpred.py", line 53, in main
    with open(sys.argv[3], 'r') as refpdbfile:
FileNotFoundError: [Errno 2] No such file or directory: 'best_qdope.pdb'
No valid current contacts.

Can you please let me know how this can be solved?

the comand line is:

run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

Thanks!

The target variable, which is set to the first part of the sequence file name, should be sequence according to your given command line. However the error message thinks it is c600m2-tail-before-lysin-sequence. Have you edited the command line for simplicity?

I'm a little worried about the floating point errors, which are almost certainly from CNS as there are 3 CNS calls and the error shows up 3 times. I would test your CNS installation to see if everything is okay - maybe try some examples from their website. An error there would explain why no output is generated.

You can ignore any errors below line 5 of the above program output for now, as they are follow-on errors. In particular the question title isn't the 'real' error here, the real error is that no PDB files are generated on line 160.

By the way, I wouldn't use 10 500 as parameters for iteration number and number of models. It will take ages to run. I'd just use the defaults (omit the arguments), which we find gives good results.

Finally, I'll point out that the DMPfold repository has an issue tracker where you can raise any future issues. Though we may as well finish debugging this one here since we've started!


Visualize network (graph) using "line" layout I want to visualise a network in a way that the nodes are along a line (and the edges are not on that line).
I have included a simple example below.

Which network visualisation programs (Cytoscape, Gephi, R igraph etc) do have this layout option?

That plot is called an Arc diagram
you can google that and use the package or program that best suits you.

Here are 2 examples of the use of Arc diagrams in R:
Example 1:

Example 2:


How I loop over this code I have 139 samples; In R for each sample I run below codes

for example for sample 1

sample_1 = whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = 1,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')

a1=(sample_1[[1]])

> class(sample_1)
[1] "list"

> a1=(sample_1[[1]])

> head(a1)
  Signature.1A Signature.1B Signature.2 Signature.3 Signature.4 Signature.5 Signature.6 Signature.7 Signature.8
1    0.5370653            0           0           0           0           0           0           0           0
  Signature.9 Signature.10 Signature.11 Signature.12 Signature.13 Signature.14 Signature.15 Signature.16 Signature.17
1           0            0            0            0            0            0            0            0    0.3383947
  Signature.18 Signature.19 Signature.20 Signature.21 Signature.R1 Signature.R2 Signature.R3 Signature.U1 Signature.U2
1    0.1123409            0            0            0            0            0            0            0            0
> dim(a1)
[1]  1 27

and for sample 2

sample_2 = whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = 2,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')

a2=(sample_2[[1]])

And finally binding a1,a2,...a139 by

sig=rbind(a1,a2)

> print(sig)
  Signature.1A Signature.1B Signature.2 Signature.3 Signature.4 Signature.5 Signature.6 Signature.7 Signature.8
1    0.5370653            0           0           0           0           0           0           0           0
2    0.7162292            0           0           0           0           0           0           0           0
  Signature.9 Signature.10 Signature.11 Signature.12 Signature.13 Signature.14 Signature.15 Signature.16 Signature.17
1           0            0            0            0            0            0            0            0    0.3383947
2           0            0            0            0            0            0            0            0    0.2837708
  Signature.18 Signature.19 Signature.20 Signature.21 Signature.R1 Signature.R2 Signature.R3 Signature.U1 Signature.U2
1    0.1123409            0            0            0            0            0            0            0            0
2    0.0000000            0            0            0            0            0            0            0            0
>

sigs.input is a common R object here

Doing these manually is too painful; Can you please help me how I do that automatically in R? I am not sure but something like looping or whatever you guess possible

Thanks

list.of.results <- lapply(seq(INSERT_TOTAL_NUMER_OF_SAMPLES_HERE),function(sample.id){
whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = sample.id,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')[[1]]
})

sig <- do.call(rbind,list.of.results)


Adding treatment groups via metadata to Seurat object? I'd like to add metadata to 6 individual Seurat objects so that after I merge the objects into one, I can later label or split by using these identifiers. For example, I'd like to append an age group and then interval across these 6 objects. For adding an interval, I've tried using the below:

AddMetaData(control, metadata = 1hr, col.name = Interval)

However, this error will occur:
Error in .AddMetaData(object = object, metadata = metadata, col.name = col.name) :
  object '1hr' not found

Appreciate any help, thank you!

You should have fed "1h" to your function call but you fed 1h instead. The latter is looking for an object with that name and as the error shows, there is no such object.


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

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


Double-counting coverage of overlapped read pairs EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region.

I used samtools depth to calculate the depth of coverage for samples in the whole Exome region using a GRCh37_ref.bed. These samples are sorted and duplicate marked. I ran this calculation on a few hundred samples to determine how much more sequencing needed to be done and found something interesting. Some samples had 100X coverage at 96,000,000 mapped reads and some samples had 100X coverage at 83,000,000 mapped reads. The samples that had 100X coverage at 83,000,000 reads had read pairs overlapping certain regions of the bedfile (read 1 was covering the same coordinates as read2 to some extent).

I then found out that samtools depth double counts these overlapped regions even though they are technically from the same molecule in sequencing and would be a duplicate bases in the read. See the image for an example:

Each base in this Duplicate Overlap region is double counted when it shouldn't be. How do you guys deal with this duplicate overlap region when calculating coverage across a bam file?

samtools (v1.7):

samtools depth -a -b exome_targets_GRCh37.bed sample1.bam > sample1_depth.txt
awk '{s+=$3}END{print s/NR}' sample1_depth.txt

samtools depth is a simplified version of samtools mpileup, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use samtools mpileup instead:

$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...

I haven't shown the entirety of the output and you'd want to use the -aa flag to ensure that all bases are returned.

Note that samtools mpileup is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to -Q 1 (the default is -Q 13, which you'd want to change).

As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use bamPEFragmentSize from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.


Can I get the graph generated by cellranger I’ve run the cellranger analysis pipeline on single cell RNASeq datasets. I can import the matrix and graph-based clusters into R. Doing this I can optimise the dimension reduction and plot cells with coloured by clusters generated by cellranger. I wish to optimise the graph-based clustering as well.

Is it possible to either:

obtain the graph generated by cellranger from the files generated. (I can only find the final clusters in the analysis output)
reproduce generating a similar graph (used for clustering) from the output matrix (preferably with R or Python.

With either approach, I am looking for a list of edges for nearest neighbors in the gene-barcode UMI matrix to use as input for a clustering algorithm.

Cell Ranger

You can't download the tSNE coordinates for cells directly from the Analysis tab of the fancy, polished .html document that Cell Ranger produces. If you have access to the machine on which the pipeline was run, you can grab the intermediate files. Everything is well documented by 10X. They explain defaults etc. on that page. The tSNE coordinates are available at:

analysis/tsne/2_components/projection.csv

and the contents look like:

$ head -5 analysis/tsne/2_components/projection.csv
Barcode,TSNE-1,TSNE-2
AAACATACAACGAA-1,-13.5494,1.4674
AAACATACTACGCA-1,-2.7325,-10.6347
AAACCGTGTCTCGC-1,12.9590,-1.6369
AAACGCACAACCAC-1,-9.3585,-6.7300

Seurat

And the second point: Cell Ranger is replicating the Seurat workflow. Seurat uses the Louvain method for community detection by default. Check out the intro tutorial and see if you can follow along. This tutorial works with 10X data, in fact the first function call in the tutorial reads in the 10X matrices. If you are already plotting in R you'll speed through this tutorial. Seurat is great to work with, you'll have tSNEs and clusters in no time. And there's so much more you can do from there!


Behavior of `--reference` flag with samtools sort I'm about to run samtools sort (version 1.6) and I'm reviewing the available configuration options.

$ samtools sort
Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read name
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    Write final output to FILE rather than standard output
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
               Specify a single output file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]

I'm curious about the --reference flag. I've never used it before with samtools sort and I'm wondering what behavior it invokes. The online documentation doesn't even mention the flag. Under what circumstances would I want/need to specify the reference genome?

You would only need to supply a reference when using CRAM output encoding. The CRAM format allows storing entries without explicitly encoding the base calls, since one can instead encode the differences to the reference genome.

With BAM output the --reference option has no effect.


Understanding the fasta File Format I'm a computer scientist teaching algorithms development in the Fall. One of the algorithms we teach is called Edit Distance, and our folklore is that it is used to compare RNA sequences (is this actually true in practice?).
I would like to have students implement the edit distance algorithm and run it on actual SARS-COV-2 sequences, so I'm trying to understand exactly what I get from the GenBank database. I downloaded this one: https://www.ncbi.nlm.nih.gov/nuccore/1798174254
I am looking at the genomic.fna file. So this is apparently the FASTA file format, and the lines that begin with >MN988669.1 ... are comments. I see comments like:
>MN988669.1 Severe acute respiratory syndrome coronavirus 2 isolate 2019-nCoV WHU02, complete genome

Followed by an RNA string. Is this the start of a new sequence for a different coronavirus specimen? So I could have students extract each of these and run edit distance and then produce a dendrogram or something? How do I find more information on where the samples came from? Is this the right file to be using, or should I use the gbff file? And are the PDB files interesting to me at all (I actually do know what PDB files are)?
Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)? Are there other useful algorithms / data that might be interesting for students to study in this vein? Specifically interesting to me would be graph-search algorithms, minimum spanning trees, and network flow. Also any NP-complete algorithms that we could run backtracking on. Obviously taking the theoretical study of algorithms to something as currently topical as the coronavirus has pedagogical value.
Thanks
EDIT:
Based on comments below here is what is taking shape.

Have students implement vanilla EditDistance (which there seems to be some disagreement about which algorithms are named what, so let's say insertions and deletions only, which I would call Longest Common Subsequence LCS). Then a variant that also does alignment (i.e. full Levenshtein distance computation, which I would call EditDistance, but Wikipedia calls the Needleman-Wunsch algorithm with a gap penalty of 1), then maybe Needleman-Wunsch with different gap penalties (if someone tells me what would make sense biologically).
Have students implement basic heirarchical clustering / phylogenetic tree generation algorithms a la https://www.ncss.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Hierarchical_Clustering-Dendrograms.pdf.
Have students run their sequence alignment variants and different clustering algorithms on SARS-COV-2 sequences and report on how the choices of parameters in 1 and 2 change the result and therefore potentially the analysis.
Ask some open ended written response questions on what this might mean for society, whether this introduces ethical considerations for an algorithms designer or are they doing just the math, etc.

My learning objectives (as they are now shaping up) are:

Students will understand that just because their algorithm comes with a proof of correctness, it doesn't mean it's the correct algorithm for the job.
Students will understand that different models / parameters to a model result in different outcomes, and thus even computational problems are not purely computational.
Doing theoretical computer science / mathematics is not devoid of ethical considerations.

I would very much appreciate thoughts on the above.

Your understanding of FASTA format is about right. The type of basic problem you're eluding to we term "sequence alignment"- edit distance might be okay for teaching but in practise we use other algorithms, e.g. you might be interested in the Needleman–Wunsch or Smith–Waterman algorithms. Richard Durbin et al. wrote a great book that covers these any much more https://www.amazon.com/Biological-Sequence-Analysis-Probabilistic-Proteins/dp/0521629713
Tracking mutations etc. requires more than just alignment though, see "phylogenetics" (i.e. building genetic trees) and "variant calling". Also checkout what the nextstrain team is doing https://nextstrain.org/ncov/global
In general, looking for practical applications for your algorithms is great, but be very careful before drawing any real world conclusions about the coronavirus outbreak from such analyses


How to maximize fastq parsing with FastqGeneralIterator (Bio.SeqIO.QualityIO) I'm contributing to a python-based project that uses Biopython to analyze fastq files. It currently uses SeqIO.parse, which populates various structures with all of the fastq information (including converting quality scores). There is apparently a faster (lighter-weight) parser called FastqGeneralIterator that doesn't populate all of these items.

I'd like to use FastqGeneralIterator, but be able to perform some of the functions SeqIO.parse offers on the fly (not for every sequence). For example, I'd like to convert base quality to PHRED scores for specific sequences, but I don't see that function available.

Is there an easy way to take a tuple output by FastqGeneralIterator and convert it into a proper Bio.SeqIO SeqRecord?

I had not come across FastqGeneralIterator before, but I will start using it for some of my work!

One answer is to replicate the code in the FastqPhredIterator, which does exactly what you are looking to do.

However, My personal preference would be to use the biopython infrastructure as much as possible, especially when it comes to the tricky aspects of record parsing.  If the records you want extra info on are uncommon, then using FastqPhredIterator for phred quality extraction might work for you:

import io
from Bio.SeqIO.QualityIO import FastqGeneralIterator, FastqPhredIterator

with open("input.fastq", "rU") as handle:
    for title, sequence, quality in FastqGeneralIterator(handle):
        # do something with title, sequence, and quality
        # Set special_case to True when you want more info on the record
        special_case = False
        if special_case:
            record_stream = io.StringIO("\n".join([title, sequence, "+", quality]))
            record = next(FastqPhredIterator(record_stream))
            # do something special with record

What we are doing here is combining the record strings into one and dumping them into a StringIO object, which implements the file handle duck-type. We then hand the in-memory handle to FastqPhredIterator to parse the FASTQ record again. Make sure to set special_case to True only when you want extra info on a record, as the special_case branch makes at least two extra copies of the FASTQ record.


Greedy Motif Search Using Probability Matrices In my greedy search with pseudocounts algorithm in my bioinformatics course, I did not follow the pseudocode since I wanted to solve the problem in my own way. Unfortunately, although my answer is correct, the answer we are supposed to give is the first correct answer, not the last.

The problem works as such: Given t strings, you make a profile matrix from the first k-mer in the first string, and then you look at the next string and using that profile matrix, you find the best k-mer match, and then you combine that k-mer with the one from the first string to make a new profile matrix, etc. You eventually do this for all k-mer windows in your first string of DNA. The way I did it was by making a probability profile matrix from the best k-mer groups in each case, and then multiplying the greatest probabilities in each column. This gives me a score, and if the next score is higher (thus, more probable) that group of k-mers is used. Below is my code:

 def highest_probability(string, k, matrix):
    score = 1
    temp_string = ""
    c = 0
    best_string = ""
    temp = 0
    while c+k < len(string):
        for i in range(c, c+k):
            if string[i] == 'A':
                score *= matrix[i-c]['A']
                temp_string += 'A'
            elif string[i] == 'C':
                score *= matrix[i-c]['C']
                temp_string += 'C'
            elif string[i] == 'G':
                score *= matrix[i-c]['G']
                temp_string += 'G'
            elif string[i] == 'T':
                score *= matrix[i-c]['T']
                temp_string += 'T'
            if i == c+k-1:
                if score >= temp:
                    temp = score
                    best_string = temp_string
        temp_string=""
        score=1
        c+=1
    print(best_string)
    return best_string

def GreedyMotifSearch(DNA, k, t):
    best_motifs = [i[0:k] for i in DNA]
    score=0
    c=0
    while c+k<len(DNA[0]):
        k_mer = DNA[0][c:c+k]
        motif=[k_mer]
        matrix = make_matrix(k_mer, matrix=[], i=0)
        for i in range(1, t):
            compare_string = DNA[i]
            good_motif = highest_probability(compare_string, k, matrix)
            motif.append(good_motif)
            matrix = make_matrix(good_motif, matrix, i)
        temp_score=1
        for dictionary in matrix:
            temp_score *= max(dictionary.values())
        print(temp_score)
        if temp_score > score:
            score = temp_score
            best_motifs = motif
        c+=1
    return best_motifs

def make_matrix(string, matrix, i):
    if len(matrix)==0:
        for ch in string:
            new_dict = {'A':0, 'C':0, 'G':0, 'T':0}
            new_dict[ch] = 1
            matrix += [new_dict]
        print(matrix)
        return matrix
    elif len(matrix)==3:
        for j in range(len(string)):
            matrix[j][string[j]] += 1
        Sum = [sum(matrix[0].values()),sum(matrix[1].values()),sum(matrix[2].values())]
        for columns in range(len(matrix)):
            for keys in matrix[columns]:
                matrix[columns][keys] /= Sum[columns]
        print(matrix)
        return matrix

DNA=['GGCGTTCAGGCA',
'AAGAATCAGTCA',
'CAAGGAGTTCGC',
'CACGTCAATCAC',
'CAATAATATTCG']
k=3
t=5
print(GreedyMotifSearch(DNA, k, t))

The output gives me ['TTC', 'GTC', 'TTC', 'GTC', 'TTC'] which is equally correct to the answer given in the problem output which is TTC,ATC,TTC,ATC,TTC. Unfortunately, my answer is the last correct occurrence. I tried adjusting my highest_probability function with if score is greater than temp instead of greater or equal, but that gives me the answer ['GGC', '', '', '', ''] which I don't understand at all. My question is 1) Why is this the case? and 2) How can I fix this? And does anyone have a simpler way of doing this?

Assuming you're talking about Implement GreedyMotifSearch from Rosalind, while I haven't actually finished this particular one, I've got an idea why TTC,ATC,TTC,ATC,TTC doesn't work. Here's my explanation from the "Questions" section of that problem:

  The profile is created from the first kmers in each string, and the most likely string from those is CAG. The answers will generally be similar to this.

I've been trying to do Rosalind by only using the information provided in the questions; this is one of the many problems where there wasn't enough information in the Question for me to arrive at a correct solution when I tried it. It's possible that the debug datasets were added to this question since I last looked at it, as they seem to describe this situation:

TEST DATASET 1:

Input:

3 4
GCCCAA
GGCCTG
AACCTA
TTCCTT

Output:

GCC
GCC
AAC
TTC

  This dataset checks that your code always picks the first-occurring Profile-most Probable k-mer in a given sequence of Dna. In the first sequence (“GCCCAA”), “GCC” and “CCA” are both Profile-most Probable k-mers. However, you must return “GCC” since it occurs earlier than “CCA”. Thus, if the first sequence of your output is “CCA”, this test case fails your code.


Is there a publicly available tumor-normal sample? I am looking for a publicly available matched tumor-normal sample. I need Illumina fastq reads (or an aligned bam file, since I could extract the reads from it) from a tumor and a matching, non-tumor control set from the same patient. Ideally, I would also have a truth set of known variants in the sample, but I can live without that.

I checked the Next Generation Cancer Knowledge Network but although they do seem to have appropriate bam files, which I could convert to fastq, all of the ones I found were "controlled access" and getting access seems very complicated.

Is there something like the public NA12878 sample for tumor-normal analyses?

Please take a look at

An open access pilot freely sharing cancer genomic data from participants in Texas

Although I have not worked on this paper but author claimed that both tumor and normal data is publicly available. In case if you still need something more, than you can work on virtual normal data. And there are several papers which provide machnine learning techniques to overcome the deficiency of normal data (not much efficient but still they give considerable results). Hope this information will help you. Thanks


Codon usage analysis for whole genomes I am new to bioinformatics. So if these questions seem you to a bit childish please forgive me.

I have two queries.

I am intending to perform a codon usage analysis followed by correspondence analyses for multiple microbial whole genomes of one bacterial species to find the association with the isolation source and ST type. Is this the right approach?
I have concatenated all the CDS in a single genome and joined all the genomes (average genes per genome is 4100). The problem is all the CDS when concatenated also have the stop codons in it. CodonW, famous for this analyses, cannot begin the analyses in the presence of stop codons. How can I remove stop codons?

Kindly suggest me a solution to the problem and concept.

This python script will remove all stop codons from your fasta file. I called it remove_stops.py

#!/usr/bin/env python3
import sys
from Bio import SeqIO
stop_codons = ["TAG", "TGA", "TAA", "UAG", "UGA", "UAA"]

fasta_it = SeqIO.parse(open(sys.argv[1]), 'fasta')
for fasta in fasta_it:
    name, sequence = fasta.description, str(fasta.seq)
    try:
        assert(len(sequence) % 3 == 0)
    except:
        raise IOError("length of sequence {} was not divisible by three".format(name))
    print(">{}".format(name))
    print("".join([sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3)
                   if sequence[i:i+3] not in stop_codons]))

You can run it like so:

python remove_stops.py stop_codons.fasta > filtered.fasta

This program won't work if there are any fasta records that do not have a length divisible by three. This indicates a sliding translation stop and these sequences should be trimmed individually. Also, all the genes must be in the 5'->3' direction.

It works by looking at each codon [sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3), then just makes sure they are not a stop codon: if sequence[i:i+3] not in stop_codons]


Kruskal-Wallis rank sum test output confusion I am running a multiple comparison using the non-parametric Kruskal Wallis test (K-W), using the ggpubr library and I am a bit confused about the results.

When i just run the KW-test using "base R" the result is different and I am not sure if there is an issue with the data or I am doing something wrong.

In case of ggpubr im getting significant result, i mean the p value is much lower than the critical threshold, but in case of base R this does not occer. I am posting the code and results figure.

Any suggestion would be highly appreciated why results are varying ..
 Data

library(ggpubr)
df1 = read.csv("NEW_RBP/CORPLOT/RBP_DISEASE_CLUSTER/RBP_DISEASE_C1_C4_AVG.txt",header = TRUE,sep = "\t")
head(df1)
my_comparisons <- list( c("HSC", "LSC"), c("LSC", "Blast"), c("HSC", "Blast") )

ex <- melt(df1, id.vars=c("gene"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

kruskal.test(HSC ~ LSC, data = df1)

    Kruskal-Wallis rank sum test

data:  HSC by LSC
Kruskal-Wallis chi-squared = 416, df = 416, p-value = 0.4908

Your problem comes from using an incorrect test, the formula parameter says:

  a formula of the form response ~ group where response gives the data values and group a vector or factor of the corresponding groups.

When you apply the kruskal  test, you don't have groups in the LSC variables! You aren't comparing the mean of several groups based on LSC!
When you use stat_compare_means it is doing a wilcox.test (it hints to it in the help page "a list of additional arguments used for the test method. For example one might use method.args = list(alternative = "greater") for wilcoxon test.").

When you use the wilcox.test we get the same p-value:

wilcox.test(df1$HSC, df1$LSC)

    Wilcoxon rank sum test with continuity correction

data:  df1$HSC and df1$LSC
W = 110980, p-value = 4.849e-12
alternative hypothesis: true location shift is not equal to 0

In addition, probably you should read about the paired argument, if this two values where measured on the same samples you should use it (but I can't tell it for sure because I don't know what does it mean each variable name).


What is and how to detect a dephased read In the methods of this paper, the authors say:

  To simplify analysis, we first removed any dephased reads in our library (last 6 bases of read did not match the expected sequence).

I have read this post and this one, and I think that a dephased read is one (according to the first paper) where the last 6 bases are not detected with good quality. In particular, I think this may be an example of a dephased read:

@...
CACACACACCAGTAAAAAANAGTTGCNNCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+
AAAAAEEAEEEEEEEEEE6#EAEEEE##E#####################################

where the quality of the read drops in the middle. Is that correct? If yes, should they be all of bad quality or at least one? And, if not, how can I detect a dephased read?

Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be.  They must be expecting some constant 6 bases at the end of each read, and they are throwing out all reads that don't have those 6 bases.  This would include N's, but a small deletion or insertion would also cause the last 6 bases to not match.  This is in theory detectable and fixable; the dropSeq protocol for instance catches instances where improperly short cell barcodes shift the UMIs, but these authors look like they decided not to try and bother, and just throw out everything that didn't match their exact expectations.


What does the yellow mean in this image from Virus Pathogen Resource? What does the yellow mean in this image? It's from https://www.viprbrc.org/brc/

Sulfur atoms are shown in yellow.

The molecular viewer that you use is JSMol (JMol ported to the web).

Atoms are colored by element: grey C, blue N, red O and yellow S.
If you wonder how other atoms would be colored, see JMol's Default element colors, by periodic table.


c++ program to find length string an char* I write under code but I can't run it when I run by g++ . erro is about  strlen.
errors:

a.cpp: In function ‘bool checker(char*, int)’:
a.cpp:185:26: error: ‘str’ is not a member of ‘std’
 185 |  for(int i = 0; i < std::str.length(word); i++)
  |                          ^~~
a.cpp: In function ‘int status(char*, int, int)’:
a.cpp:192:40: error: ‘str’ is not a member of ‘std’
   192 |  for(int i = 0; i < 4 && nw*4+i < std::str.length(pmt[pn]); i++)
  |                                        ^~~
a.cpp:199:10: error: ‘str’ is not a member of ‘std’
  199 |  if(std::str.length(pmt[pn]) > nw*4)
  |          ^~~
a.cpp: In function ‘int main()’:
a.cpp:225:22: error: ‘str’ is not a member of ‘std’
  225 |     len_words = std::str.length(pmt[i])/2;
    |                      ^~~

I think you're confusing C++ and C strings, which are very different. When you have a char *, that's a C-style "string". It's really just a bunch of bytes with a zero byte at the end, demarcating the end of the string. To find the length of a C string you include <cstring> (or <string.h> in C) and call strlen, like this:

/** Print every character in the string */
void print_chars_in_string(char* word) {
    for (size_t i = 0; i < strlen(word); ++i) {
        putchar(word[i]);
    }
}

strlen is actually pretty simple though, you can write your own version as an academic exercise like this:

size_t deluxe_bespoke_strlen(char* str) {
    size_t length = 0;

    while (*str++)
        ++length;

    return length;
}

Since the while-loop condition will be true until you hit the null terminator (the zero byte at the end), this is how you would implement the strlen function without a for-loop, otherwise you would need the strlen function for the strlen function.

In C++, strings are actual objects, with methods you can call. To get the length of a C++ string simply do this:

std::string my_name = "Jose";
std::cout << my_name.length() << std::endl;

If you have a C string, you can convert it to a C++ string with the built-in constructor:

const char* my_name = "Jose";
std::string my_name_cppstring = std::string(my_name);

You can also get the C version of a C++ string by calling the c_str() member function like this:

std::string sequence = "GATTACA";
size_t length = strlen(sequence.c_str());

My recommendation is to not mix C and C++. C is my favorite programming language, but that doesn't mean everyone should write everything in C. Since you're using C++ I'm guessing that's your preference, so I would suggest using C++ strings instead of C strings, just remember to include <string>.

As for your specific question, the compiler is giving you the correct message: str is not a member of std. You don't add the namespace to variables. This is how that works:

#include <cstdlib>
#include <iostream>

namespace CustomDeluxeBespokeStandardNamespace
{
    size_t custom_bespoke_strlen(char* str)
    {
        size_t length = 0;

        while (*str++)
            ++length;

        return length;
    }
}

int main(int argc, char *argv[])
{
    if (argc == 1) {
        std::cerr << "No input(s)." << std::endl;
        std::exit(EXIT_FAILURE);
    }

    while (*++argv) {
        std::cout << *argv << ": " << CustomDeluxeBespokeStandardNamespace::custom_bespoke_strlen(*argv) << std::endl;
    }
}

So if you run it like this:

testprogram horse house mouse cat

You should get this output:

horse: 5
house: 5
mouse: 5
cat: 3

So to recap, you were adding the namespace to your variables, and the compiler was giving you an error because it correctly could not find anything by that name in the std namespace. Even if it had though, C style strings do not have member functions. There are no such things as objects in C, at least not in the sense that you can define member functions for an object, so that was the second error.

As an aside, you will note all my examples using strlen both custom and actual return size_t. This is not an accident; in both C and C++ the return type of a call to strlen is size_t. You have to be careful because in your code you compare an int to a size_t. Comparing a signed value to an unsigned value is Not Good because it will work only if the signed value can be faithfully represented after a conversion to an unsigned value. If it can't though, you're in for a bad time.


LD analysis in PLINK based on reference and a SNP list I'm new to PLINK and genetics, and getting confused with two PLINK commands for LD analysis:

plink --bfile hapmap --r2 --ld-window-r2 --ld-snp-list --ld-window
plink --noweb --bfile --clump hapmap --clump inputfile --clump-field P --clump-p1 --clump-p2 --clump-r2 --clump-kb

It seems to me that for #1, the input files are a reference data and a SNP list file, while for #2, it is a reference data and a .assoc file.

In the output for #2, p-values are included with other results (according to PLINK manual) which is not the case with the output from #1. If I'm interested in finding proxy SNPs for the SNPS on my list in hapmap data, will it be sufficient for me to use the output from #1?

After invoking command #1 using PLINK 1.9, its output gave me 6 columns: "Chromosome number", "chromosomal position", "SNP name" for query/target SNP and 1 column of "R2". I cannot carry out command #1 because I do not have the .assoc file and only have a list of SNPs.

My main question is - what is the difference between these two commands and which of the two should give me the best result for finding proxy SNPs for my list of SNPs.

It looks like the first command is calculating pairwise linkage disequilibrium r2 between your SNPs.

The second command is clumping your dataset based on a p-value over a file of association test results, which I guess is a pre-cursor to computing polygenic risk scores.

Clumping is the process of thinning your SNP list to only keep uncorrelated SNPs so that each representative region has only one marker. In this case you are asking PLINK to keep the marker with the highest p-value.

I'm not sure what your definition of proxy SNPs is but if you are looking for tag SNPs to a SNP of your interest, you should go with command 1 and choose a proper LD r2 threshold.


Interpreting large z scores from microarray data I have microarray z scores for gene set enriched (GSE) pathways that range approximately from -3 to 14. I understand that an absolute value of 1.96 has a p-value = .05 and significant for our purposes. My question is whether I can consider extreme z scores as more significant?

This is a statistical question. What a higher z value means is that it is more extreme (if you assume a normal distribution), thus a more extreme value is less likely to happen by chance. Which is translated to a (lower) p-value.

What it is significant and not is arbitrary. Prior to the experiment you must decide which is your threshold, how much chance you allow. Apparently you decided to select 0.05 as your threshold, which means that one of 20 results might be due to chance. But your extreme z score is not more significant, it is less expected (under the null hypothesis, those p-values might change if you change your null hypothesis, ie you think your genes expression don't follow a normal distribution) but as significant as your 0.0499999 p-values.


umap highlighting two different models I'm trying to create a umap for single cell data from human samples and ptx samples. I can get the umap to where it shows the umap with the different clusters but I want to show where the ptx samples are and where the human samples are.

My code is as follows:

#create the Seurat object

OD_10K_HUMAN <- CreateSeuratObject(counts = HUMAN_OD_10K.data, min.cells = 0, project = "human")
SD_5K_HUMAN <- CreateSeuratObject(counts = HUMAN_SD_5K.data, min.cells = 0, project = "human")
BNL.5K <- CreateSeuratObject(counts = SD_BNL_5K.data, min.cells = 0, project = "ptx")
BNL.10K <- CreateSeuratObject(counts = OD_BNL_10K.data, min.cells = 0, project = "ptx")
BNM.10K <- CreateSeuratObject(counts = OD_BNM_10K.data, min.cells = 0, project = "ptx")
BNM.5K <- CreateSeuratObject(counts = SD_BNM_5K.data, min.cells = 0, project = "ptx")

#merge data
scData <- merge(BNL.10K, y = c(BNL.5K, BNM.10K, BNM.5K, SD_5K_HUMAN, OD_10K_HUMAN), add.cell.ids = c("A", "B", "C", "D", "E", "F"), project = "HTB2876")

mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = scData), value = TRUE)
length(mito.genes)
scData[["percent.mt"]] <- PercentageFeatureSet(scData, pattern = "^MT-")

scData[["log_nCount_RNA"]] <- log2(scData[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
scData <- subset(scData, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
scData <- NormalizeData(scData, normalization.method = "LogNormalize", scale.factor = 10000)
scData <- FindVariableFeatures(scData, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(scData)
scData <- ScaleData(scData, features = VariableFeatures(object = scData), vars.to.regress = c("nCount_RNA"))
scData <- RunPCA(scData, features = VariableFeatures(object = scData))
#DimPlot(scData, reduction = "pca")

numPC = 20
scData <- FindNeighbors(scData, dims = 1:numPC)
scData <- FindClusters(scData, resolution = 0.4)
scData <- RunUMAP(scData, dims = 1:numPC)

DimPlot(scData, reduction = "umap", label = TRUE)

By default, DimPlot colors the cells by their "identity class", stored in Ident(scData) in your case. This indeed contains the clustering information.

However, you can you use the group.by argument for DimPlot to pass another column from scData@meta.data to color cells by. In your case you would need a metadata column that stores the ptx/human information.

The merge function automatically stores the original "project" information of the individual Seurat objects that you combine under orig.ident in the metadata. Since you set the project names to either "ptx" or "human", the orig.ident column should contain the information you need.

Something like this should work:

DimPlot(scData, reduction = "umap", group.by = "orig.ident", label = TRUE)


Q: What analyses can I perform with a completely phased genome assembly? Phasing a genome is the process of determining which variants lie on which copy of each chromosome. For example, position 5,430,500 might have an A on the paternal chromosome, but a T base on the maternal chromosome.

A haplotype block is a section of the genome for which we know which haplotype has what variants. A completely phased genome has phase blocks that are the full length of each chromosome. Every variant on every chromosome can be assigned specifically to the maternal to or the paternal haplotype.

My question is, what can we learn from having completely phased genomes that we cannot from an unphased genome, or from a phased genome with short haplotype blocks?

Some population genetic inference methods require the data to be phased because they look for long stretches of matching alleles between individuals, rather than just matching genotypes. Here, 'genotype' data is un-phased, and 'allele' data is phased.

This is because there is a degree of ambiguity when it comes to matching heterozygous positions - if the data is un-phased, you don't know if 2 individuals share mutations on the same chromosome (here, I mean different copies of the same chromosome - like either copy of chromosome 1, for example). This is important, because the fact that the mutation is shared on the same chromosome means that they are much more likely to share recent common ancestry with one another.

Let's say we have genotype data for two individuals (where the number corresponds to the number of non-reference alleles present):

Ind 1

0 ----- 1 ----- 2 ----- 0

Ind 2

0 ----- 1 ----- 2 ----- 0

If were looking at unphased genotype data only, we know the alleles match exactly in the 1st, 3rd and 4th positions, because they are homozgyous (either 0 or 2 copies of the non-referencee allele), but because the second position is heterozygous, we can't be sure which chromosome the non-reference allele is on. This is because the haplotype configurations for individual 1 could be either

0 ----- 0 ----- 1 ----- 0 Chromosome1a
0 ----- 1 ----- 1 ----- 0 Chromosome1b

or

0 ----- 1 ----- 1 ----- 0 Chromosome1a
0 ----- 0 ----- 1 ----- 0 Chromosome1b

i.e., we don't know which chromosome the non-reference allele at position 2 is on.

If we then phased the data, and found that the non-reference allele was on chromosome 1a for Individual 1, and chromosome 1b for Individual 2 there is less evidence of recent shared ancestry than otherwise. I think this is more or less the approach that 23andme takes when they are inferring your relationship to other individuals in the dataset. After writing this, I realised there is a much better description of what I just said on the 23andme website!

For example, the above applies to the ChromoPainter/fineSTRUCTURE algorithms. This means that they are able to detect more subtle genetic differences between individuals than if you were to use non-phased data.

This is a slightly niche example, and there occasions where you need phased data, but this is the one I am most familiar with :). Off the top of my head, methods such as Multiply Sequential Markovian Coalescent (MSMC), RELATE, and Extended Haplotype Homozygosity all require the data to be phased.


How to get rid of the temp files while using "featureCounts" for extracting readcounts from bam files? I used this command for extracting read counts for almost 60 samples using featureCounts.

featureCounts -a /gencode.v27.primary_assembly.annotation_nochr.gtf -F GTF -p -s 2 -T 8 -o counts.txt sample.sorted.bam

I see many "temp-core-100373-A4BF0116E148.sam-TH00-BK004414.tmp" files like this dumped in my directory. I'm able to delete them only manually. Is there a way to get rid of these tmp files from my directory while using the command?

You shouldn't delete them while featureCounts is running, it's creating read-name-sorted temporary files so it can more easily look at both mates together. Just wait for it to finish running and it will clean up the files itself. Doing so beforehand is liable to cause an error.


europmc and tidypmc R library for extracting or making metadata from publication I'm trying to use this europmc r libray where I have a list of pmids to look for. I tried with pubtator but its bit complicated.In Europmc i can all the annotated terms etc.
library("europepmc")

example list of PMIDS
30024784
30555165
30510081
31688884
31516032
28588019
29286103

Now what I'm doing is looking each ID using epmc_details function which is not the way i would do if i have to look for hundreds.
epmc_details(ext_id = '30510081')

I have to question how can I run the epmc_details in a loop or some other way where i can look for PMIDS one by one and save the result in a data frame.
The epmc_details is returned as a list. The structure of the list is as such
[1] "basic"           "author_details"  "journal_info"    "ftx"             "chemical"        "mesh_topic"      "mesh_qualifiers"
[8] "comments"        "grants"

I would only like to save basic,chemical,mesh_topic,mesh_qualifiers in one data frame.
For example if my first id is this 30510081 the dataframe should have first column as my ID which is basically basic[1]and rest of the information appended to the next columns.
such as
 ID chemical mesh_qualifiers mesh_topic gene

Any suggestion or help would be really and highly appreciated
I was looking at europmc site through the browser
this was one of my query
when i highlight the keyterms i do see in the abstract itself all the keyterms are getting annotated but when i do the same query search through R I do see empty results as such why there is a difference?
$chemical
# A tibble: 0 x 0

$mesh_topic
# A tibble: 0 x 0

$mesh_qualifiers
# A tibble: 0 x 0

I found better way of getting data from pubmed using the tidypmc library.
library(tidypmc)
doc <- pmc_xml("PMC6365492")
doc
txt <- pmc_text(doc)
txt

count(txt, "section")

cap1 <- pmc_caption(doc)
filter(cap1, sentence == 1)

tab1 <- pmc_table(doc)

sapply(tab1, nrow)

tab1[[1]]
attributes(tab1[[2]])

collapse_rows(tab1, na.string="-")

library(tibble)
x <- xml_name(xml_find_all(doc, "//*"))
tibble(tag=x) %>% count("tag")

library(tidytext)
x1 <- unnest_tokens(txt, word, text) %>%
  anti_join(stop_words) %>%
  filter(!word %in% 1:100)
#  Joining, by = "word"
#filter(x1, str_detect(section, "Case description"))

filter(x1, str_detect(section, "Results"))
count(a$word)

tbls <- pmc_table(doc)
map_int(tbls, nrow)

tbls[[1]]

collapse_rows(tbls, na.string="-")

But if i understand it can use one PMC id at a time. Again keeping my original question how can i put this in a loop to query lets say i have 100 PMCID and get it result and store in a dataframe.
After using tidypmc i found that i parse all the publication based on attributes or tags. Such as title,abstract, results etc etc.
Lets say I'm interested in the table tags information where they have metadata of patients as well as others.So if a paper contain multiple tables I would like to store each of them in a data-frame under the respective publication.Since I have multiple IDs to search and do the save as above mentioned. How to put this through a loop or can it be done without using loop ?
Any suggestion or help would be really appreciated as always.

The idea for this code is to first convert PIDs to PMCIDs, then run the tidypmc in a loop over the PMCIDs. The only problem is that tidypmc failed to retrieve tables from most of the IDs in your example list.
library(tidyverse)
library(tidypmc)
library(httr)
library(jsonlite)
example_pids <- c(30024784, 30555165, 30510081, 31688884, 31516032, 28588019, 29286103) %>% as.character()

#-- Convert to PMC ids
convertPIDtoPMCID <- function(pids) {
    #-------- Make API request
    pids4query <- paste(pids, collapse = "%0D%0A")
    idconv_req <- paste0("https://www.ncbi.nlm.nih.gov/pmc/utils/idconv/v1.0/?ids=", pids4query, "&idtype=pmid&format=json&versions=no&showaiid=no&tool=&email=&.submit=Submit")
    pids_json <- GET(idconv_req) %>% content("text") %>% fromJSON()
    #-------- Get info from JSON
    pid2pmc <- pids_json$records %>%
        select(pmcid, pmid) %>%
        as.data.frame()
    rownames(pid2pmc) <- pid2pmc$pmid
    pmcids <- pmcids <- pid2pmc[pids, "pmcid"]
    return(pmcids)
}
example_pmcids <- convertPIDtoPMCID(example_pids)

#-- Try to get data with tidypmc
pub_tables <- lapply(example_pmcids, function(pmc_id) {
    message("-- Trying ", pmc_id, "...")
    doc <- tryCatch(pmc_xml(pmc_id),
                    error = function(e) {
                        message("------ Failed to recover PMCID")
                        return(NULL)
                        })
    if(!is.null(doc)) {
        #-- If succeed, try to get table
        tables <- pmc_table(doc)
        if(!is.null(tables)) {
            #-- If succeed, try to get table name
            table_caps <- pmc_caption(doc) %>%
                filter(tag == "table")
            names(tables) <- paste(table_caps$label, table_caps$text, sep = " - ")
        }
        return(tables)
    } else {
        #-- If fail, return NA
        return(NA)
    }
})
names(pub_tables) <- example_pids

#-- Inspect results
pub_tables$`30555165`$`Table 1 - Patient demographic and baseline characteristics`
pub_tables$`29286103`$`Table I - Sample summary.`

Tables will require quite a bit of tidying after this. Good luck!


Is outlier detection the right choice for determining if a value lies singificantly away from the mean in a population in gene frequency? I have the following dataframe that is comprised of enzyme substrates as the variable and then the unique number of times enzymes capiable of degrading this substrates appeared in my genomes and then tot total count of times all of those enzymes appeard in my samples. For example, a 7 unique enzymes were found to be capiable of degrading chitin, and they occured a total of 298 times across my genomes.

I wish to see if any particular substrate has significantly larger occurance in both unique and total, than the rest of the substrates. Would outlier detection be appropriate for this or would a simple T-test do?

I agree you have the right idea because that data doesn't look at all linear, so regression analysis looks really complicated. Your approach a straight comparison keeps it simple and is within the biological question you are asking,

However, I would use a non-parametric such as Wilcoxon instead of a T-test. I don't fully know the null distribution you are comparing it against. Wilcoxan doesn't require an assumption of normality. T-tests invariably require a transformation to normality and that gets complicated, but a log or sqrt transformation might be okay. If Wilcoxan produces a no significant difference you report it and move onto a T-test.

It is difficult to use outlier analysis with 9(?) observations, thats used if you have a way-out result and can afford to exclude it (>3 standard deviations). Outlier analysis is to remove a single result that is messing up your analysis and doesn't affect the bulk of the data. You alot of observations to be considering outlier analysis.


RNASeq: Normalization, stabilization, gene length and rlog I was thinking about the best method for normalization, which takes gene length into account (in order to compare genes)...

Do you think I can do that? :
- taking raw counts and dividing each gene by its length
- using the function rlog (DESeq2) on these counts divided by gene length (I would modify the rlog function to allow it to be used on decimal data).

I was wondering if it would be ok to do that.

Thank you very much for your help

I am not sure what you mean by "decimal data", however, if you really want to use "gene length" information for normalization, take a look at EDASeq. Unlike DESeq or EdgeR, its normalization step takes gene length into account to produce "normalized counts", which you can feed into DESeq or EdgeR according to their manual.

Below is from the said manual:

  The normalized counts (or the original counts and the offset) obtained
  using the EDASeq package can be supplied to packages such as edgeR
  (Robinson, McCarthy, and Smyth 2010) or DESeq (Anders and Huber 2010)
  to find differentially expressed genes.


Masking sites in a vcf file I need to mask all sites in a vcf file flagged by the 1000 Genomes Project as being unfit for population genetic analyses. The sites for all chromosomes are available at:
1000Genomes masked sites
From above mentioned link you can get both bed files and Fasta files.
I wanted to mask all sites using VCFtools but VCFtools is not longer being updated and therefore, can't me used on the my vcf file because you can only VCFtools upto VCF file format v.2. And file format for my file is v.3. Can someone suggest another tool that I can use. Insights will be appreciated.

Have a look at BCFtools, which I think does most of the stuff that VCFtools does. I notice there's a --mask option in BCFtools consensus, which may do what you want:
-m, --mask FILE
 BED file or TAB file with regions to be replaced with N
 (the default) or as specified by the next --mask-with
 option. See discussion of --regions-file in Common
 Options for file format details


Performing t.test to obtain p.values for each gene I have FPKM values for Normal(80 cases), Group1(135 cases), Group2(147 cases), and Group3(102 cases) for 160 genes (Groups were obtained by clustering tumor data).

I want to perform t-test (two-sided, unpaired, unequal variance) between Normal and Group1, Normal and Group2, and Normal and Group3 such that I get p.value for each gene.
The t.test function in R returns single p.value and not gene-wise (I don't want to do t-test in Excel).

Which package/function should I use to do this in R?

You can simply apply() t.test() to your matrix. In general, though, I expect lmFit() from the limma package will be simpler to use.

Here is an example with t.test():

m = matrix(runif(100, max=c(rep(1, 50), rep(1.5, 50))), nrow=10) # make some data
apply(m, 1, function(x) t.test(x[1:5], x[6:10])$p.value)

Neither t.test() nor lmFit() should be expected to produce the same results as excel, since excel is inferior in every way and should never be relied upon for any statistical analysis (plus lmFit is doing something much much more complicated than a simple t-test).

As an aside, FPKMs are highly problematic for use in statistics. If possible, I strongly encourage you to get either per-gene counts or per-gene estimated counts and to use those in limma or another differential expression package instead, as the results will be much more reliable. I presume, though, that you're using TCGA data, so you may be stuck with FPKMs.


Targeted hybrid capture bam files showing reads outside target regions in my bed file Sorry if this is a stupid question I am new to analysing targeted hybrid capture data.
I am analysing my HTS data from a targeted hybrid capture enrichment method and subsequent sequencing on a NextSeq.
Looking at my bam files, I see there are reads that are overlapping into my non target regions.
For example - my .bed file has two target regions in chr17:

However I still have reads overlapping between 7578739-7578918 (following example is position 7578850):

Quite a lot of these reads don't hit any part of the targeted regions in my bed file. Are these mistakes? How have they been pulled down for sequencing when they aren't in the target regions of the baits?

It strongly depends on exactly what your protocol is. Note that there are a good number of tools for QCing capture data (here, here).
Off-target reads are a standard observation in capture experiments; recall that it is target enrichment, not target perfect purification, as commenter swbarnes2 suggests.
The real question is whether the proportion of those off-target reads is too high for your application, or otherwise suggests problems. Those depend on details of the technique that you will likely have to figure out yourself.
Here are some tables from the first of those QC papers (supplement), showing on-target proportions for various capture samples (sorry, it would have been a bit of effort to get text of the tables rather than images):

These data suggest that on-target reads are 60%-80% of mapped reads. Note that they also lose a lot of reads in mapping, so we cannot exclude the possibility of mis-mapped garbage reads contributing to off-target mapping.


How to extract exome on-target reads in batch? I was given a list of target regions in BED and many exome alignments in BAM. I was asked to extract on-target alignments from these BAMs to save disk space. I know I can use bedtools to extract sub-BAMs. I am thinking to write a script to apply bedtools to all BAMs at my hand, but I speculate there may be some more convenient command lines to achieve my goal. How would you do? Thanks.

It sounds like bedtools intersect will work for you:

bedtools intersect -wa -a <alignment.bam> -b <region.bed> > <intersect_alignment.bam>

  When a BAM file is used for the A file, the alignment is retained if overlaps exist, and exlcuded if an overlap cannot be found.  If multiple overlaps exist, they are not reported, as we are only testing for one or more overlaps.

It's easy enough to wrap this up into a for loop for batch processing (line breaks can be removed if desired):

 for alnBam in alignmentFiles*.bam;
   do echo ${alnBam};
   bedtools intersect -wa -a ${alnBam} -b <region.bed> > intersect_${alnBam};
 done

There are a lot of other options which should cover most related uses. For more details, just run "bedtools intersect" (with no arguments).


How I deal with this expression set? I have a gene expression raw counts like below for 16 patients

         A  B    C      D   E   F
A2M     511 1623 665    208 553 469
AADAT   34  137 372     7   52  124
ABCB1   119 114 123    22   25  89
ABCB11  27  186 200    27   30  49
ABCC2   1   10  21     1    8   3

For each patients I also have some characteristics

Patients                          IC          TCs
A1                                3            3
B1                                2            0
C1                                2            0
D1                                0            0

I want to find genes differential expressed in these patients regarding these characteristics but I don't know how to do that

Any help please?

Look up a tutorial on DESeq2.  To keep things simple for your first time, I'd reduce your phenotypes to "high" and "low" bins instead of trying to use the cell count numbers as they are.  Then once you understand what the software is doing, you can get fancier.


Getting this information for GISTIC2 Long time I am struggling how to provide GISTIC2 required input like Num_Probes and Segment_Mean

I have called somatic copy number by scatngs and I have copy number output file (attached);

https://www.dropbox.com/s/r6nef2v1j2di6q6/my.copynumber.caveman.csv?dl=0

This is few lines of scatngs copy number output

1,1,13116,665401,2,1,3,1
2,1,668374,96324627,2,1,2,1
3,1,96326095,96419709,2,1,1,0

In scatngs reference paper says in this file we have

Segment number

Chromosome

Start position (origin-1)

End position (origin-1)

Major copy number—normal

Minor copy number—normal

Major copy number—tumor

Minor copy number—tumor

scatngs also returns another output like below

https://www.dropbox.com/s/7i28vooixdp06c2/WTSI-OESO_121_1pre.copynumber.txt?dl=0

Chromosome  Position    Log R   segmented LogR  BAF segmented BAF   Copy number Minor allele    Raw copy number
rs62635286  1   13116   -0.0290741822638953 0.245246748556562   0.25    NA  3   1   2.35806652809599
rs75454623  1   14930   0.402869048014209   0.245246748556562   0.3095  NA  3   1   4.14031506947036

I am wondering how I can get Num_Probes and Segment_Mean From this files?

People say

Segment_Mean = log2(tumour copy number / 2)

So, does that mean I should add up Major copy number—tumor and

Minor copy number—tumor to get tumour copy number? for getting Num_Probes what should I do?

EDITED

  people says

  scatngs uses read counts at SNPs to generate calls and so the
  locations of the SNPs is what I need to
  generate the number of markers in each segment. (eg each SNP = 1
  marker).

  Do they mean, for instance if the position of SNP  rs62635286 is
  13116, so the num_prob is 13116?

Any help please?

Thanks a lot in advance

Depends on the software to generate the segmentation file;

For example I have used scatngs for segmentation that for sample returns

Segment number
Chromosome
Start position (origin-1)
End position (origin-1)
Major copy number—normal
Minor copy number—normal
Major copy number—tumor
Minor copy number—tumor

By  log2(Major copy number—tumor +  Minor copy number—tumor /ploidy of each sample) you would get segment mean

For number of probes

The number of probes depends on where the data came from.
1) From SNP arrays:  Number of array probes within each segment.
2) From exon capture data:   Probes are the ends of exons and so the number of probes are equal to the number of exons  X 2.
3) From whole genome sequencing:  Probes are the boundaries of bins used to generate the copy number data and so the number of probes are equal to the number of bins +

For example in my case for scatngs

Each SNP is a marker, and so you have to count how many of these SNPs are on each segment

I have asked this question in

Extracting the number SNP in each range


Identifying nested/overlapping reading frames in nucleotide sequences I'll start by saying I'm a total stranger to this field, so I'd love to be corrected if I misunderstood something along the way.
I'm writing a python code to identify all the reading frames in a sequence of nucleotides. To my understanding, a reading frame is defined by a start codon and a stop codon, which made total sense to me until I learned you can have nested/overlapping reading frames.
My question is:
Assuming the start codon "AUG", and the stop codon "UAA":
AUG AAA GGG AUG CCC UUU UAA AGA CUC UAA
How do you identify the reading frames? Is each stop codon referring to a single start codon? If so, how can you tell to which one?
I tried looking this up and could only come up with articles explaining the effects of overlapping reading frames, rather than how are they interpreted. I would greatly appreciate even a nudge in the right direction.
Please let me know if I was not clear enough or if I got something wrong.
Thanks!
EDIT: After reading the answers here and researching further I learned that ORF identification on DNA sequences isn't as simple as finding a single start and stop codon, and is complicated beyond the scope of my question. I'll probably default to using online tools to do this from now on.

An RNA sequence has six possible reading frames depending on the offset (0, 1, 2) and the direction (forward, reverse complement). All frames can potentially be used for  protein translation, but only one frame per protein (i.e. there is no shifting going on when translation is happening, excluding shifts introduced by splicing).
If you need to find all the possible protein sequences within an RNA sequence, there are plenty of programs already available (e.g. getorf and transseq from the EMBOSS tool kit).


Read and write FASTA files with more information than id and sequence I am trying to read a fasta file, manipulate is in Python (using BioPython) and then write it back.

The format of my sequences is like:

>k119_5 flag=0 multi=141.0706 len=473
AGGTTAGTCAGCACCGTTTCCGTGGTGCTGCCTTTCGCTTCAAAACCGACGGCGTCTATTACTGCATCCACGCCCCGGTGACCTGCCGTTTGTTCAATAATTGACTGTGCCGGATCGCTGTCTTCATCAAAATTAATCGGGATCGCGCCGTAGCGGTCGGCGGCGAAATGCAAGCGGTAGGGATGATGATCAACAACAAAAATCTGTTCCGCACCGAGCAACCGTGCACAGGCGATTGTCAACAATCCCACAGGACCAGCACCATAGACTGC
AACGCTTGAACCTTGTTGGATCTGCGCATTTTTTGCTGCCTGCCATGCCGTTGGCAGAATATCAGAAAGGAAAAGCGCTTTATCATCTGAAAGCAAAGGCGGTACTTTAAACGGCCCCACATTCCCTTTAGGGACGCGGACATATTCCGCCTGCCCACCAGGAACGCCGCCATACAGGTGACTATAACCAAACAATGCCGC

which is the result of the assmebly using Megahit. While trying to read, it only recognizes k119_5 as the id, while I want the whole first line as the id of the sequence. Is vital information I want to keep!

Any ideas how to work with this?

The description field in the SeqRecord object has the information you are looking for:

>> from Bio import SeqIO
>> s = SeqIO.read('genome.fasta', 'fasta') # single sequence fasta file
>> s.id
k119_5
>> s.description
k119_5 flag=0 multi=141.0706 len=473

Edit: as an aside, if you write a SeqRecord object using SeqIO's write method, the FASTA header should contain the whole original header.


Chromosomal order of supplementary alignment in BAM file "SA" tag I have several long-read sequences, which when aligned to a human reference genome, aligns to multiple chromosomal fragment as a result of chromothripsis shattering the linear DNA and randomly piecing them together. This creates for a single primary alignment multiple supplementary alignments.

With a single example of a read (only showing parts of all the bamfile information) being

cdaf1fcb-27bf-4890-8ac7-a097dfa12c38    16      chr11   15151487        60

With the primary alignment starting at chr11 at 15151487

Then the supplementary alignment can be seen as follows
SA:Z:chr13,49827969,+,1839S956M25D1584S,60,53;chr11,15152459,-,833S749M22D2797S,60,72;chrX,138324656,-,2540S404M2D1435S,60,20;chr11,15152066,+,278M11D4101S,60,20;

Where the alignments are to chr13, 11, X and 11. With both regions on chr11 being close to the initial region from the primary alignment, making the regions on chromosome 11 the same (I know I haven't calculated the full chromosomal interval using the Cigar string for this example).

I assume the order of that supplementary information above is in the actual order. But the problem is that my DNA should be circular which doesn't fit with the order of the alignments since the circle would then be chr 11 (primary, start of circle) -> chr 13 -> chr11 (same as primary) -> chrX -> chr11 (back to the start of circle)

My question is:
Q) So is the order of the supplementary information just pieces together randomly or is it true that I can assume the order is the true order of the sequencing read?

So thanks a lot for your help.

There is nothing in the SAMtags specification that dictates the order of the entries in an SA tag.  It may be entirely dependent on the aligner that you used.

Looking at my own data that was aligned with bwa, I cannot see a definite pattern.  Longer matches tend to be first, but not always.

Your aligner documentation may have more information on how it uses SA tags.


How to obtain the standard gene name I have a list of genes and proteins extracted from literature. How can I convert them into standard gene names.

Gene name

Using an Ensembl ID, Uniprot ID or NCBI ID it is much easier to look up your gene. So I would strongly reccomend using one of those as a key. Personally, I'd use Uniprot, but Ensembl is probably better if your genes are human.
One reason for this is also that different resources will use different names, especially for humans.

There are a few considerations, so I will first say how this could be done properly...

Custom Uniprot matcher with weights

Some names are used by different protein, so which one is which would depend on context, co-occurrence with another gene.

Your dataset might be extremely noisy. Therefore, you likely have assigned a confidence score to your matches, yes? It would be nice to give a score to the name based on whether it is the first name of the protein or the fourth synonym in the list.

Therefore to do a conversion properly you really ought to make your own conversion system based upon context and Uniprot data. Uniprot contains the most data on synonyms and notes and binding partners etc. The dataset/API is in XML and Biopython does not have a parser, but I have some notes on the matter and can share some code if needed as I have done something similar (e.g. JSON for human gene/protein names to Uniprot IDs).

Lazy match

The NCBI entrez API is rather flexible and good at presenting to you matches for a strange names. So I would suggest using that.


Official documentation for textual BLAST output? This is pretty vexing: I got asked to explain the meaning of the match line in a blastp output, especially the subtle difference between a blank and plus. E.g.:

 Query  181  NHIEEL...
             N I+EL...
 Sbjct  181  NQIKEL...

For me this is under the "everyone knows this" category, but I wanted to have an authoritative reference I could cite.

Alas, none of the documents on the BLAST help page seems to have this information. My Google fu also failed as I was not able to find an authoritative source which I could cite.

Also, I'm not sure whether BLAST introduced this notation or it simply picked it up somewhere else.

Any pointers appreciated.

I think gringer is correct to say this is because of the BLOSUM62 matrix and a ‘+’ symbol inidcating a positive score for a mismatch but here is a little more explanation.

BLASTP looks for (nearly) identical short protein matches pairwise, called seeds, and then extends the seeds in each direction one residue at a time using arbitrary gap penalties for indels and a matrix for substitutions. Based on some thresholds and more, the highest scoring alignments are kept and ultimately reported.

The default choice of matrix for BLASTP is BLOSUM62, which looks like this. Obviously matches have high positive scores, but also mismatches have varying positive or negative scores.

If you look up your changes:

H -> Q; ; score 0

E -> K; +; score 1

Note that this matrix is symmetric so the scores for H -> Q == Q -> H and E -> K == K -> E.

In terms of why E -> K has a positive score, well not all amino acid changes are viewed equally. BLOSUM matrices are not based on any biological theory but instead use the log-odds ratio, a likelihood approach, that estimates from real sequence alignments how frequently different amino acid changes occur. The more often E -> K substitutions are observed in these alignments, the higher the substitution will score in the derived matrix.

Hopefully these scores are intupretable based on the biochemical properties of the amino acids. For example, R -> K scores 2, a  high score for a mismatch, but both amino acids are basic and so share a similar structure, presumably making a change more likely to occur and/or more likely to maintain the protein function and so not be harmful.

For your case, I don't have a good biological explanation for why E -> K scores positive, another author notes this is surprising, but ultimately the reason it is given a positive score is just because regardless of why, E -> K does apparently frequently occur, at least in the alignment data used to derive this scoring matrix. Incidentally, BLOSUM matrices are pretty ancient and based on few alignments by modern standards, but nonetheless appear to work quite well.