text:
Bioinformatics flowchart for Materials and Methods
I am going to start work on the Materials and Methods for a paper I contributed to a lot in the bioinformatics processing part. I have seen some papers nowadays tend to first present the processing steps that they followed as a flow chart, to show the reader what happened generally before going into details in the written M&M section.
Are there any specific guidelines for creating this flowchart? Do you have any tips or advice from personal experience?

(I am aware that this may be off-topic. I have a meta question about it here)

Still waiting for the PLOS Computational biology "10 simple rules..." on this one ;)

I don't know... Is a flow chart actually better than a concise description possibly divided into sections with headers and bullet points? The flow chart you linked doesn't strike me as particularly clear or informative given it takes a full page with a big figure and a wall of text. (I'm not being critical I'm genuinely asking...)

If you're using a workflow management system for pipelining like snakemake or ruffus, you can actually build a flowchart based on the rules/functions in your code:
Snakemake example with the --dag flag:

Ruffus example using pipeline_printout_graph:

You'll probably still need edit the image aftwards in software like inkscape or illustrator though.

I think it's better and easier to understand to make a more concise flowchart than this, and to just publish the code along with the paper. The first flowchart seems to list every single input on which  the analysis was run. The second one is far too complex to read.

@Pallie agreed, the flowcharts can get quite complicated fast, or contain elements that you don't really need to show. You can edit them as suggested, but it would be nice to be able to display just certain parts of the pipeline in the flow chart as an option.

Flow charts in vogue, when they are not ... reviewers are very quick to point this out


text:
How to retrieve dataset names from python statsmodels?
The Python statsmodels library grants access to many R-datasets which can be downloaded e.g. with:
import statsmodels.api as sm
df = sm.datasets.get_rdataset(dataname='iris', package='datasets').data

But for this you have to know the id or the name of the dataset.
sm.datasets.get_rdataset??
Signature: sm.datasets.get_rdataset(dataname, package='datasets', cache=False)
Source:
def get_rdataset(dataname, package="datasets", cache=False):
    """download and return R dataset

    Parameters
    ----------
    dataname : str
        The name of the dataset you want to download

There is a list of available datasets here, but this list does not contain the identifiers. Also not the detailed view.
Is there an overview of all the available datasets that contains the ids e.g. "iris"?
There is an index that is supposed to list all the datasets, but e.g. 'rice' as a name does not work.

Here is the website of the creator of rdatasets I guess with the info you are looking for:
http://vincentarelbundock.github.io/Rdatasets/datasets.html
Here is a short look at it:
Package     Item    Title   Rows    Cols    n_binary    n_character     n_factor    n_logical   n_numeric   CSV     Doc
boot    acme    Monthly Excess Returns  60  3   0   1   0   0   2   CSV     DOC
boot    aids    Delay in AIDS Reporting in England and Wales    570     6   1   0   0   0   6   CSV     DOC
boot    aircondit   Failures of Air-conditioning Equipment  12  1   0   0   0   0   1   CSV     DOC
boot    aircondit7  Failures of Air-conditioning Equipment  24  1   0   0   0   0   1   CSV     DOC
boot    amis    Car Speeding and Warning Signs  8437    4   1   0   0   0   4   CSV     DOC
boot    aml     Remission Times for Acute Myelogenous Leukaemia     23  3   2   0   0   0   3   CSV     DOC
...

I assume the Item column are the identifiers.

As stated above, that does not seem to work. Try for example 'rice'.

Did you change the package parameter to DAAG for the dataset rice? Since it’s not in the package datasets according to the index site.

You can view all the available datasets at this URL https://github.com/vincentarelbundock/Rdatasets/blob/master/datasets.csv
The available ones in statsmodels are those with the tag datasets which looks like line 372 - 455.
I didn't test them all, but everyone I tried from a small random sample seemed to work fine.

Please, try this code to get all modules in datasets:
import statsmodels.api as sm
import inspect

def get_list_of_modules_in_datasets():
  # get list of strings for the objects in datasets
  list_str_in_datasets = sm.datasets.__all__
  print(list_str_in_datasets)
  # get list of objects in datasets
  list_objs_in_datasets = [eval('sm.datasets.' + dsstr) for dsstr in list_str_in_datasets]
  print(list_objs_in_datasets)
  # get list of modules in datasets
  list_modules_in_datasets = []
  for dso in list_objs_in_datasets:
    if (inspect.ismodule(dso)):
      list_modules_in_datasets.append(dso)
  print(list_modules_in_datasets)
  print("Number of Modules in datasets = %d" %(len(list_modules_in_datasets)))
  return (list_modules_in_datasets)

if __name__ == '__main__':

  list_modules_in_datasets = get_list_of_modules_in_datasets()

  # print data for all dataframes from the datasets
  for ds in list_modules_in_datasets:
    df = ds.load_pandas().data
    print(list(df.columns))
    print(df.shape)
    print(df)

source: Stack Exchange

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

sort: votes


text:
Java error when trying to run UBCG jar (java.lang.ClassNotFoundException)
I wanna install UBCG.JAR in ubuntu 18.04 but I cannot run the java File UBCG.JAR.
when I type the code java -jar UBCG.JAR an error message appears.
Exception in thread "main" java.lang.NoClassDefFoundError: org/json/simple/parser/ParseException
at java.base/java.lang.Class.forName0(Native Method)
at java.base/java.lang.Class.forName(Class.java:398)
at org.eclipse.jdt.internal.jarinjarloader.JarRsrcLoader.main(JarRsrcLoader.java:56)

Caused by: java.lang.ClassNotFoundException: org.json.simple.parser.ParseException
at java.base/java.net.URLClassLoader.findClass(URLClassLoader.java:471)
at java.base/java.lang.ClassLoader.loadClass(ClassLoader.java:588)
at java.base/java.lang.ClassLoader.loadClass(ClassLoader.java:521)
... 3 more
$ java -version
Openjdk version "11.0.7" 2020-04-14

Thank you very much the problem solved, it was about JDK version

Could you please either accept an answer, or write up (and accept) your own solution to this problem?

You are probably using the wrong java version it should be 8
try "java -version"

I solved this problem but UBCG.JAR does not work with java 11 . It is necessary to install openjdk-8-jre and switch the default java.
Install openjdk-8-jre
$ sudo apt-cache policy openjdk-8-jre
$ sudo apt install openjdk-8-jre=1.8.0_292

Switch the default java
$ sudo update-alternatives --config java
*            /usr/lib/jvm/java-8-openjdk-amd64/jre/bin/java   1081

Checking the default java
    $ java -version

openjdk version "1.8.0_292"
OpenJDK Runtime Environment (build 1.8.0_292-8u292-b10-0ubuntu1~18.04-b10)
OpenJDK 64-Bit Server VM (build 25.292-b10, mixed mode)
Run the UBGC.jar
$ java -jar UBCG.jar
      ###############################
      #         UBCG ver. 3.0       #
      ###############################

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


text:
Finding simple sequence from reads with significant overlap
I wrote a "script" to pull out reads from a huge fastq file in an iterative manner, by finding homology to the previous sequence.
It should be relatively easy to overlap them and assemble the gene, however traditional assemblers like SPAdes fail due to obviously low coverage.
I am hesitant to just feed it more faked reads, when I know I have all the data necessary to just overlap them into a single sequence.
The reads are set up like this:
GATCAAACATC
     AACATCAGTTAG
         TCAGTTAG
              TAGAGGATAGC

What kind of assembler is out there where I can feed it this kind of terribly-formatted data?

Miniasm will probably work for this; it's designed to be an ultra-fast, greedy assembler, which has been used for people doing real-time assemblies from MinION data. Because of its greedy nature, it can misassemble regions where the overlapping fragments are ambiguous, so it's worth using some other method in addition to miniasm to confirm that the assembly is correct.
By default it looks for overlap regions that have a coverage of 3 reads when joining (option -c), but I assume that could be reduced down to 2 reads, as in your example (which has a minimum overlapping coverage of 2).

Thanks for the suggestion. This actually wound up working in a roundabout way!


text:
bash loop for strelka2 processing for multiple files from one directory
In the below bash loop for strelka2, 3 bam files are being stored in a variable. The two questions are:
Is this proper way to loop through multi bam and wait for workflow to be created for strelka2 processing,
and can you create a workflow in a directory without the bam in it. Thank you :).
contents of dir
s1.bam
s2.bam
s3.bam

bash
dir=/path/to/bam
data=/path/to/processed/data

cd $dir
 for bam in ${dir}/*.bam; do sample_basename=$(basename "${bam}"); done # loop through dir and get bam

${path_to_strelka}/bin/configureStrelkaGermlineWorkflow.py \
  --bam ${dir}/${sample_basename} \
  --referenceFasta <fasta> \
  --runDir ${data}
  wait  # wait until workflow in data is created before executing below  in data
${data}/runWorkflow.py -m local -j 20

The loop you're using does not work, as it overwrites the value of sample_basename each time. Plus, given that you're passing in the full path to the bam (using $dir/$sample_basename) anyway, I don't see the reason to use the sample_basename variable at all.
Just this:
dir=/path/to/bam
data=/path/to/processed/data
ref_fa=/path/to/ref/fasta

cd $dir

for bam_file in *.bam
do
  ${path_to_strelka}/bin/configureStrelkaGermlineWorkflow.py \
    --bam ${bam_file} \
    --referenceFasta ${ref_fa} \
    --runDir ${data}
  ${data}/runWorkflow.py -m local -j 20 #I'm assuming this goes inside the loop. If not, move it outside the loop. You don't need a `wait` in a shell script that's working line-by-line already.
done

Thank you very much :).


text:
protein secondary structure matching/comparison scoring
I want to compare two secondary structures of aligned proteins. I do not know exactly how to do it well.
Example:
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC

EEEHHHHHHHTTTTTTTTCTTHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHEEEECCCCCCCCCCEEEEEEEEEEEEEE

Example of matching:
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC
|||        |||||| |||                    |||||||||     ||||||      |   |||    ||
EEEHHHHHHHTTTTTTTTCTTHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH----HHEEEECCCCCCCCCCEEEEEEEEEE

Hypothetically: 50% aligned - > score 0.5 - >  thats scoring is wrong:
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC
||||||||||||||||||||||||||||||||||||||||||
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHEEEEEEEEEEEEEEEHHHHHHHHHHHHHHHHHHHHHHHH

This case will have similar score...
This is a puzzle for me, I don't know how to approach it technically.

Normally 3D structures of protein are compared and the RSMD is reported. How did you get your SS, with JPred or from PDBs? Comparing SS might result in some biases: if one of the two protein is getting used as a template for the other say or if both are from PDBs but the SS was called different (if so, use DSS on both). Also note that JPred is like 75% accurate at best.

@MatteoFerla I suspect that these are predicted secondary structures.

@James Exactly my worry (I see a helix of a single AA in a coil). The solution is straightforward, but I am not sure how one could add error estimates to said solutions, which seems like a must...

I suppose using data from JPred, PSIPRED, SPIDER2 and doing the 9 pairwise alignments would suffice on a push.

@MatteoFerla I agree. In lieu of a structure, usually, a meta-prediction would suffice. I think MTG should add why the sequences are being aligned. I think different tools should be used depending on if they are checking to see how much secondary structure prediction error is in a sequence compared to a reference sequence, or if they are checking secondary structure as a proxy for structural divergence.

It is a predicted SS with GOR method. I try to add sorting by SS to my algorihtm. First version will use garnier (emboss) method. Sequences are aligned because program align similar sequences. Then I add SS to check aligned proteins (additional validation). I plan to use spider or jpred in future.

In summary, this is a very provisional analysis, providing a starting point to further de novo prediction methods. We all assume the crystal structures are not present here

Hey i solve this simply by use of SOV alghoritm. https://www.researchgate.net/publication/319526145_A_modified_definition_of_Sov_a_segment-based_measure_for_protein_secondary_structure_prediction_assessment
cheers


text:
Gene size and power of detection in sequence divergence analyses
I want to look at violation of the strict molecular clock and compare dn/ds ratio to test differences in evolutionary rates in a set of genes. Some of these genes are really short (30AA), is it likely to prevent the detection of rate changes ?
Thanks

If your genes are really short you may want to find a mathematical hack to get a confidence interval (several published), otherwise the jitteriness will be too much.

I can't find anything about size correction/confidence interval in this type of analysis. As dn/ds is a ratio, why would it be biaised with shorter genes? I feel this can be an issue but I can't find anything explianing why and how to deal with that.

Would concatenating several genes make the trick? For example, based on the evolutionary model that fit them the best (concatenate all genes under GTR+I model for example) or all genes related to the same function?

It is not a bias (no size correction needed), it is just that a ratio of two integers gets 'jittery' with low values (say 1/1 = 1, 1/2 = 0.5!). Concatenating different genes would assume the selection pressures are the same, so that would be not advisable. If I am completely misunderstanding and you were going to average (or similarly process) these anyway, averaging them will get rid of the jitters. In terms of errors the simplest and most rudimentary calculation is a margin of error, but that is pointless if the values are to be combined.

Ok thanks. I'm not willing to combined the values if it's avoidable.


text:
Density plot, scale it to 0-1
Doubt regarding density plot what is the scale being plotted in the Y axis.Is it possible to scale it to 1?.
This is my code
w <- read.csv("Normal_Myeloid_Dev_stages/Myeloid_non_coding_non_CDS_CORRVALUE.txt",header = TRUE)
head(w)
#names(w)[1] = "Sample"
head(w)

#w.plot <- melt(w)

df <- w

df_melt=melt(df,id.vars="Sample")
#head(df_melt)
tail(df_melt)
df_melt$Group <- gsub('[0-9]', '',df_melt$variable)
head(df_melt)

p1 <- ggplot(aes(x=value, colour=Group), data=df_melt)
p1 + geom_density()

ggplot(df_melt, aes(x = value,y=..density.., fill = Group)) + geom_density(alpha = 0.5)+
  theme_bw(base_size=30)+
  theme(axis.text.x=element_text(angle = 45, size=45, face="bold", hjust = 1),
        axis.text.y=element_text(angle=0, size=50, face="bold", vjust=0.5),
        plot.title = element_text(size=40, face="bold"),
        legend.title=element_blank(),
        legend.key.size=unit(1, "cm"),      #Sets overall area/size of the legend
        legend.text=element_text(size=30))

dev.off()

Second question regarding finding the overlapping area between various condition let say I have four condition and i want to find out what is the overlapping or intersecting area between them .
I tried on of the stack solution firstly is it the right way to do it for finding overlap? if yes can i try finding all the condition overlap such as I have HSC,CMP,GMP and Mono.
FDensity = approxfun(density(df_melt$value[df_melt$Group=="CMP"], from=.80, to=1))
MDensity = approxfun(density(df_melt$value[df_melt$Group=="Mono"], from=.80, to=1))
plot(FDensity, xlim=c(.80,1), ylab="Density")
curve(MDensity, add=TRUE)

FminusM = function(x) { FDensity(x) - MDensity(x) }
Intersect = uniroot(FminusM, c(.8, 1))$root
points(Intersect, FDensity(Intersect), pch=20, col="red")

integrate(MDensity, .80,Intersect)$value +
  integrate(FDensity, Intersect, 1)$value

Any suggestion or help would be really appreciated
My dataframe small subset
dput(head(df))
structure(list(Sample = structure(c(9L, 10L, 11L, 12L, 1L, 2L
), .Label = c("CMP1", "CMP2", "CMP3", "CMP4", "GMP1", "GMP2",
"GMP3", "GMP4", "HSC1", "HSC2", "HSC3", "HSC4", "Mono1", "Mono2",
"Mono3", "Mono4"), class = "factor"), HSC1 = c(1, 0.901758000194052,
0.880971458703505, 0.900098712568466, 0.873507054094022, 0.906315987662777
), HSC2 = c(0.901758000194052, 1, 0.945122894369186, 0.955453795442752,
0.844777154811663, 0.960260413763721), HSC3 = c(0.880971458703505,
0.945122894369186, 1, 0.931891528389177, 0.829521543809914, 0.93636611655036
), HSC4 = c(0.900098712568466, 0.955453795442752, 0.931891528389177,
1, 0.852155096935692, 0.947192034704188), CMP1 = c(0.873507054094022,
0.844777154811663, 0.829521543809914, 0.852155096935692, 1, 0.88406635204624
), CMP2 = c(0.906315987662777, 0.960260413763721, 0.93636611655036,
0.947192034704188, 0.88406635204624, 1), CMP3 = c(0.883085447599108,
0.928929790213059, 0.937278860333014, 0.923139729112196, 0.870595169861622,
0.958330253497026), CMP4 = c(0.903500691840647, 0.931015449611016,
0.915510499719426, 0.949501343959892, 0.89151869084162, 0.960206335031192
), GMP1 = c(0.859320159223793, 0.833454007185579, 0.821479398953591,
0.837605565046685, 0.856892209387879, 0.8612685942226), GMP2 = c(0.894680165774456,
0.944129985346879, 0.92461599550388, 0.937108867595306, 0.875930650641323,
0.966512850930959), GMP3 = c(0.881416749993221, 0.925001567707159,
0.923015420025178, 0.920758302664238, 0.863998116268766, 0.95069438729762
), GMP4 = c(0.887230702792805, 0.911614382740949, 0.898704725226052,
0.938310378760322, 0.875956580458754, 0.93865124308053), Mono1 = c(0.80834274546097,
0.803706049167148, 0.785535924609134, 0.794438848049474, 0.794402126200663,
0.810382054101301), Mono2 = c(0.829928035079747, 0.868021129791583,
0.848929312833536, 0.852953184941131, 0.818074585324121, 0.880490870699251
), Mono3 = c(0.834028389107401, 0.859818911060746, 0.852694253707916,
0.85039434449057, 0.82251614206182, 0.872740655210769), Mono4 = c(0.832524159512094,
0.858444209302146, 0.84397886802844, 0.862850406374668, 0.818436260136434,
0.870293804485003)), row.names = c(NA, 6L), class = "data.frame")

I'm not fully to have understood your second questions, but first, if you want to scale the density to 1, you can use y= ..scaled.. in your aes:
library(tidyverse)
df %>% pivot_longer(., everything(), names_to = "Variable", values_to = "Value") %>%
  ggplot(., aes(x = Value, y = ..scaled.., fill = Variable)) + geom_density(alpha = 0.4)

Now regarding the second question, if you are looking for the overlap between each of your density, you can use the package overlapping as specified in the post you are referring: https://stackoverflow.com/questions/41914257/calculate-area-of-overlapping-density-plot-by-ggplot-using-r?noredirect=1&lq=1
library(overlapping)
library(lattice)
x <- list(D1 = D1, D2 = D2, D3 = D3, D4 = D4)
overlap(x, plot = TRUE)

And if you want to extract the value of overlap, you can get it:
OV <- overlap(x)$OV

> OV
    D1-D2     D1-D3     D1-D4     D2-D3     D2-D4     D3-D4
0.5971760 0.3171644 0.7566528 0.4926057 0.6683083 0.3756224

Is it what you are looking for ? Or did I misunderstood your question ?
Dummy Dataset
# Dummy example:
D1 = rnorm(100, mean =1, sd = 2)
D2 = rnorm(100, mean = 2, sd = 2)
D3 = rnorm(100, mean = 3, sd = 1)
D4 = rnorm(100,mean = 1.5, sd = 1.5)

df = data.frame(D1,D2,D3,D4)

EDIT: Using reproducible example provided in the question
For plotting geom_density of your data:
library(tidyverse)
df %>% pivot_longer(., -Sample, names_to = "Variable", values_to = "Value") %>%
  mutate(New_Var = gsub("\\d","",Variable)) %>%
  ggplot(., aes(x = Value, y = ..scaled.., fill = New_Var))+
  geom_density(alpha = 0.4)

And to find the overlap:
d<- df %>% pivot_longer(., -Sample, names_to = "Variable", values_to = "Value") %>% mutate(New_Var = gsub("\\d","",Variable))
a <- split(d$Value, unique(d$New_Var))

g <- overlap(a, plot = TRUE)$OV

I don't think there is an option in the overlap function to scaled each overlap plot to 1 but these plot are using ggplot and so, I imagine you can extract their points and scaled them by plotting them separately in ggplot2.
Hope it answers your question

excellent you have got my question precisely let me give it a try on my data set

a little doubt for the second part since I have multiple biological replicate ,since you made a list x <- list(D1 = D1, D2 = D2, D3 = D3, D4 = D4) as you have each individual column which you are overlapping .I will attach a small part of my dataframe to show my doubt.

Since in the density plot I making a group of each replicate and plotting them as one but if i try your second part I will have to test individual sample, now my question is can i overlap it as different group rather than individual samples .

I managed to split my dataframe into list which is where i was stuck by doing this .       a <- split(df_melt$value, unique(df_melt$Group)) but im not sure if it would be an issue or not as i see the scales of my Y axis are not scaled to 1 in case of those overlapping plot

I edited my answer to show the use of your example dataset. However, I'm not sure it is possible to scaled each overlap plot to 1 (probably as you can extract individual data for each plot). Also, as you don't have massive difference in the amplitude of each density, I'm not sure the scaled function is really required for that.

excellent answer thank you for making it clear

I have a doubt what is the double peak in case of CMP GMP , I see GMP two split peaks what are those how to explain?

Could potentially be two subpopulations of the GMP group...

okay ..so here what im taking is correlation of different sample type so I have 4 replicate each in case of GMP if i understand then you are saying one of the biological replicate or replicates might have different correlation pattern ?

Yes, either one or two of your biological replicate is different from others, or you have two subpopulations in all of your replicate. I think, you can try to plot their densities separately to see how values of each replicate is distributed. Either, they are all single peak but with different means (which after merging create two peaks), or each replicate have two peaks.

that's really good suggestion I will do and verify


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

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


text:
“IndexError: list index out of range” in snakemake
I am trying to establish pipeline in snakePipes and while running DNA mapping I got next IndexError error.
(snakePipes) ~ $ DNA-mapping -i /Volumes/HD2/snakePipes_CHIP/ChiPseq -o /Volumes/HD2/snakePipes_CHIP/ChiPseq1 --mapq 5 -j 10 --dedup hg38
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 10
Job counts:
    count   jobs
    5   Bowtie2
    5   FASTQ1
    5   FASTQ2

rule origFASTQ2:
    input: /Volumes/HD2/snakePipes_CHIP/ChiPseq/SRR6761499_R2.fastq.gz
    output: originalFASTQ/SRR6761499_R2.fastq.gz
    jobid: 75
    wildcards: sample=SRR6761499

RuleException in line 9 of /Users/v/miniconda3/envs/snakePipes/lib/python3.6/site-packages/snakePipes/shared/rules/FASTQ.snakefile:
IndexError: list index out of range
Error: snakemake returned an error code of 1, so processing is incomplete!

When I am looking in FASTQ.snakefile:
rule origFASTQ1:
      input:
          indir+"/{sample}"+reads[0]+ext
      output:
          "originalFASTQ/{sample}"+reads[0]+".fastq.gz"
      shell:
          "( [ -f {output} ] || ln -s -r {input} {output} )"

line9: rule origFASTQ2:
      input:
          indir+"/{sample}"+reads[1]+ext
      output:
          "originalFASTQ/{sample}"+reads[1]+".fastq.gz"
      shell:
          "( [ -f {output} ] || ln -s -r {input} {output} )"

I do not understand why I see this error. I hope somebody can point me on the source of the error.
When I use --local parameter I do not see this error, however I get another error:
ln: illegal option -- r

usage: ln [-Ffhinsv] source_file [target_file]
       ln [-Ffhinsv] source_file ... target_dir
       link source_file target_file
    output: originalFASTQ/SRR6761495_R1.fastq.gz
    shell:
        ( [ -f originalFASTQ/SRR6761495_R1.fastq.gz ] || ln -s -r /Volumes/HD2/snakePipes_CHIP/ChiPseq/SRR6761495_R1.fastq.gz originalFASTQ/SRR6761495_R1.fastq.gz )
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

as far as I understand -r option creates symbolic links relative to link location
in /Volumes/HD2/snakePipes_CHIP/ChiPseq I have:
Bowtie2         SRR6761495_R1.fastq.gz  SRR6761495_R2.fastq.gz
command.sh      SRR6761499_R1.fastq.gz  SRR6761499_R2.fastq.gz
ChIPdict.yaml   SRR6761497_R1.fastq.gz  SRR6761497_R2.fastq.gz
deepTools_qc    SRR6761498_R1.fastq.gz  SRR6761498_R2.fastq.gz
fetch.sh        SRR6761502_R1.fastq.gz  SRR6761502_R2.fastq.gz
filtered_bam
Sambamba
rename_files.sh

I am working on macOS Catalina in bash shell and ln -r is valid option (I installed coreutils)
Thanks


text:
nanopore - where to retrieve information from the basecaller used
Where can I find information about the basecaller used in a specific nanopore run (we have a Gridion machine)?
I know already that we have Guppy v. 3.0.3, which I found out by clicking on the about tab on the MinKNOW GUI).
I was wondering, though, if additional information about the basecaller is outputted in any of the files generated when the run is finished? I checked already the final_summary.txt and it does not contain such information. I also checked the head of the fastq files and no further details.
My ultimate goal for the moment is to run medaka on a draft genome assembly. But for that I'd need to specify the medaka model. From medaka_consensus -h I see that the options available and defaults are:
-m  medaka model, (default: r941_min_high_g330).
        Available: r941_min_fast_g303, r941_min_high_g303, r941_min_high_g330, r941_prom_fast_g303, r941_prom_high_g303, r941_prom_high_g330, r10_min_high_g303, r10_min_high_g340, r941_prom_diploid_snp, r941_min_high_g340_rle.
        Alternatively a .hdf file from 'medaka train'.

Therefore I already know that my option should contain min and g303 (format of the models available described here). But know I'd need to find out the {pore} and {caller variant}

Pore? Say, "R9.4.1"? This is the pore is the physical protein (mutant curlin G with helicase) on your chip that you used. It likely will be r941, unless you are trialing out the new R10 pores.

Yes I also believe the pore would be R9.4.1. Not sure about the caller variant now

I basecall separately with guppy. In the output folder specified by --save_path or -s there are a whole bunch of .log files.
$ ls -l *.log | head
-rw-r--r-- 1 tom tom 5242714 Dec  3 11:04 guppy_basecaller_log-2019-12-02_22-02-36.log
-rw-r--r-- 1 tom tom 5242718 Dec  3 11:06 guppy_basecaller_log-2019-12-02_22-04-38.log
-rw-r--r-- 1 tom tom 5242730 Dec  3 11:08 guppy_basecaller_log-2019-12-02_22-06-41.log
-rw-r--r-- 1 tom tom 5242879 Dec  3 11:11 guppy_basecaller_log-2019-12-02_22-08-50.log
-rw-r--r-- 1 tom tom 5242811 Dec  3 11:13 guppy_basecaller_log-2019-12-02_22-11-06.log
-rw-r--r-- 1 tom tom 5242817 Dec  3 11:15 guppy_basecaller_log-2019-12-02_22-13-28.log
-rw-r--r-- 1 tom tom 5242880 Dec  3 11:18 guppy_basecaller_log-2019-12-02_22-15-47.log
-rw-r--r-- 1 tom tom 5242816 Dec  3 11:20 guppy_basecaller_log-2019-12-02_22-18-09.log
-rw-r--r-- 1 tom tom 5242777 Dec  3 11:22 guppy_basecaller_log-2019-12-02_22-20-30.log
-rw-r--r-- 1 tom tom 5242742 Dec  3 11:25 guppy_basecaller_log-2019-12-02_22-22-47.log

The information you need for medaka is at the top of the first log file (alphabetically, and by creation date).
$ head -n 13 guppy_basecaller_log-2019-12-02_22-02-36.log
2019-12-02 22:02:36.255793 [guppy/message] ONT Guppy basecalling software version 3.4.1+ad4f8b9
config file:        /opt/ont/guppy/data/dna_r9.4.1_450bps_hac.cfg
model file:         /opt/ont/guppy/data/template_r9.4.1_450bps_hac.jsn
input path:         raw/minion
save path:          test_341/minion
chunk size:         1000
chunks per runner:  512
minimum qscore:     7
records per file:   4000
num basecallers:    1
gpu device:         cuda:0
kernel path:
runners per device: 4

If you ran live basecalling and you don't have those log files you could also try a recursive grep from the output folder, e.g.
grep -nr "^model file" .


text:
Merging files gives memory allocation error
I have two files like
> head(a)
  chr position
1   1   136962
2   1   562020
3   1   672948
> class(a)
[1] "data.frame"
>

And
> dim(b)
[1] 1855235      10
> head(b[,1:3])
            SNP chr position
1:  rs62635286          1    13116
2:  rs75454623          1    14930
3:    rs806731          1    30923
4: rs200943160          1    49298
5: rs116400033          1    51479
6: rs141149254          1    54490
> class(b)
[1] "data.table" "data.frame"
>

I want to find corresponding SNP in file b matched to position in file a so I did like this that gives memory allocation error.
> unique(merge(a,b,x.by=a$position,y.by=b$position))
Error: cannot allocate vector of size 137.2 Gb

I also tried rstudio cloud but uploading the bigger file never completes
Can you help me with task please?
The expected output would be something like
rs62635286  1   13116
rs376723915 1   91515
rs147061536 1   92858
rs62642131  1   135982
rs371474651 1   158006
rs201293782 1   665401
rs138476838 1   668374
rs1401137   1   2063094
rs3128293   1   2065339

EDITED

> setDT(a)
> setkey(a, chr, position)
> setkey(b, chr, position)
Error in setkeyv(x, cols, verbose = verbose, physical = physical) :
  some columns are not in the data.table: chr,position
> setkey(b, Chromosome, Position)
> a[b, nomatch=0]
Error in bmerge(i, x, leftcols, rightcols, roll, rollends, nomatch, mult,  :
  Incompatible join types: x.chr (factor) and i.Chromosome (integer). Factor columns must join to factor or character columns.
>

You could try something like this :
First you should work with the same class of object : data.table for fast merging operation
# Turning a to a data.table
setDT(a)
# b is already a data.table

Then you might have to set the keys of both table :
# I added chr in case you need to match the same position on a different chromosome.
setkey(a, chr, position)
setkey(b, chr, position)

Finally perform the joining operation :
a[b, nomatch=0]

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

Thank you I added the output of your command in my main post. I obtained error please can you have a look?

This error means the type of your chromosome columns don't match. You can check the type of all columns with str(a). I suggest you keep the the integer format of the b table. This should update a column from factor to integer : a[, chr := as.integer(chr)]. It may work then !


text:
how to change the PC use in the dimplot and feature plot
I would like to know how to change the PC use in the dimplot and featureplot by using Seurat. for we can get the x-axis and the y-axis like PC-1 and PC-2, if I want to use PC-4 and PC-5. How to modify the code?

You can choose which PC dims to plot by specifying the dims argument within the functions. To use PC-4 and PC-5,
DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "pca")

and
FeaturePlot(object = pbmc_small, features = "Xist", dims = c(4, 5), reduction = "pca")

Seurat documentation is available here.

Hi Kohl, thank you for your advice. But after I tried that code, I got the error:"Error in Embeddings(object = object[[reduction]])[cells, dims] :
subscript out of bounds
Calls: DimPlot
Execution halted"    Do you have any idea about this?

From the docs: **reduction** Which dimensionality reduction to use. If not specified, first searches for umap, then tsne, then pca. I imagine you haven't run UMAP or tSNE yet. Try providing the reduction arg, I've updated my answer.

Hi Kohl, thank you for your advice. And I have successfully got the PC change, but I would like to know whether we can change the UMPA use? like the PC use?

If I understand your question correctly, replace reduction = "pca" with reduction = "umap". You'll also have to keep more dims when you do RunUMAP(), the default is to keep two. Be sure to add dims.keep = 5L to RunUMAP() if you want to access dims 4 and 5. If you need any more detail, you should ask another question as this has branch from your original question and would be more helpful for future seekers if it were on it's own.


text:
Running IBD using PLINK --parallel (v2)
I have merged two huge .gen files with roughly 500k samples in total - around 3 TB.
I wish to run an IBD (identity by descent) check to identify if there are identical participants in the two cohorts before doing any further analyses.
In the past I have simply used the --genome argument in PLINK and it would work fine however, in this instance the file is so large that doing such would cause the following error to appear: Error: Out of memory.  The --memory flag may be helpful.
I then opted to use PLINK v2 as it has a --parallel argument that, apparently, can be used alongside --make-king and/or --make-king-table to conduct this IBD check. I have had some trouble understanding the documentation surrounding these aforementioned arguments and would like some clarification as to how to ideally use them. Please assume that I am running on a multi-node, 24 core cluster.
Here is the current code I am running:
plink2 --gen testo.gen \
--sample testo.sample \
--make-king \
--make-king-table \
--parallel 1 24

And here is what the PLINK documentation says about the arguments for this command:
--parallel <1-based current job index> <total job pieces>
Now, I have failed to understand what the first argument does, even after playing around with the parameters a little and extensively searching the web for an answer.
What would changing this "1-based current job index" do? What does it mean? And for "total job pieces", will the adjustment of this parameter alter how many cores the jobs are distributed to or is it just an identifier as to how many compartmentalized sets are being made?
If not, is there any way to specify which nodes you would like to run a job on and specify the number of cores as well? If yes, could you please provide an example?
Thanks.

Here is an extensive documentation : https://www.cog-genomics.org/plink/1.9/parallel

Two other notes: if you’re just trying to detect duplicates, you can use —king-table-filter to greatly reduce the output file size (and you can also afford to select a small random sample of higher-MAF variants, instead of brute-forcing the entire dataset); and if you’re running on a single machine, PLINK 2.0 doesn’t need —parallel, unlike PLINK 1.9 it’ll automatically adapt to the memory limit.


text:
exclude uncultured and environmental samples from NT databases
I'm using blast locally against NT databases and I'm excluding the uncultured and environmental samples by the taxonomy ID. Is there another way to exclude them entirely, without using the taxonomy ID?

https://bioinformatics.stackexchange.com/questions/7384/taxon-exclude-list-for-searching-local-blast-database-using-blastn

If trying to gain database quality by excluding uncultured and environmental, RefSeq may be helpful.  Database can be downloaded locally. update_blastdb.pl --showall shows pre-formatted databases (refseq_genomic).  See link, "Reference Sequence (RefSeq)" at https://www.ncbi.nlm.nih.gov/guide/all/ for more detail.


text:
Given a .gb file and a locus - how to get relevant annotations in Python?
Given a .gb file and a specific locus in the genome - how can I parse the relevant annotations in Python (i.e., annotations that include that locus)?

I import the features using:
SeqIO.read('my_gb_file.gb', 'gb').features

I then wish scan/filter sequence meta data.

Is there a function or method in Biopython that does this?

Biopython is the main package for this. It's only a few lines, so it is not reinventing the wheel.
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).

A nice handy way to do operations on the features is to convert it to a pandas dataframe.
from Bio.SeqFeature import SeqFeature
from Bio import SeqIO
from Bio.Seq import Seq
import pandas as pd

        def convert(feat: SeqFeature, flatten:bool=True):
    if flatten:
        clean = lambda v: v[0] if isinstance(v, list) else v
        qualifiers = {k: clean(v) for k, v in feat.qualifiers.items()}
    else:
        qualifiers = feat.qualifiers
    return {'id': feat.id,
            'type': feat.type,
            'location_strand': feat.location.strand,
            'location_start': feat.location.start,
            'location_end': feat.location.end,
            **qualifiers}

transcript = SeqIO.read('.gb', 'genbank')
feats = pd.DataFrame([convert(feat) for feat in transcript.features])

This means that in a jupyter notebook one can display the table to look at it (feats as last line in a cell or IPython.display.display(feats)), export it to a csv (feats.to_cvs('.csv')) or subset it:
subfeats = feats.loc[(feats['type'] == 'exon') & (feats['location_start']  > 100)& (feats['location_end']  < 500)]

This is not a homework question. A list in python comes with the builtin function index. It felt to me a bit like implementing list.index, so I asked. This seems to me like an extremely common operation, so I thought it is probable that such a function already exists, and I just failed to find it. Thanks anyway for your answer!

Sorry, I am really bad at figuring out when it's a homework question and when it isn't...

No harm done. And clearly you did your best to help :)

I removed my statement about no homework and copypasted a chunk of code to convert it to a pandas dataframe which makes life easier.


text:
Comparing factor levels in deseq2
I have sample information (sampleinfo) in this format:
SampleID    condition   type
MP_1    A_1 Paired-End
MP_2    A_1 Paired-End
MP_3    A_1 Paired-End
MN_1    A_2 Paired-End
MN_2    A_2 Paired-End
MN_3    A_2 Paired-End
WP_1    B_1 Paired-End
WP_2    B_1 Paired-End
WP_3    B_1 Paired-End
WN_1    B_2 Paired-End
WN_2    B_2 Paired-End
WN_3    B_2 Paired-End

I want to compare A_1 with A_2, B_1 with B_2 ans also want to compare (A_1+A_2) with (B_1+B_2).
This is how I used DESeqDataSetFromMatrix
ddsMat <- DESeqDataSetFromMatrix(countData = seqdata[,rownames(sampleinfo)], colData =  sampleinfo, design = ~ condition)
ddsMat <- DESeq(ddsMat)

res.05_FC0<-results(ddsMat,contrast = c("condition","A_1","A_2"),alpha=0.05,lfcThreshold = 0)
res.05_FC0<-results(ddsMat,contrast = c("condition","B_1","B_2"),alpha=0.05,lfcThreshold = 0)

How can I compare (A_1+A_2) with (B_1+B_2) using the contrast argument in results function?
Of course I can make a new sampleinfo table, but I want to use the same sampleinfo as above. I tried to read about numeric contrast, but I am unable to use it correctly.

You can use contrast = c(0.5, 0.5, -0.5, -0.5), which means "the average of the first two coefficients minus the average of the second two.


text:
Obtaining HGDP project data in fasta format
I need to obtain sample data from modern humans in fasta format. I just need some megabytes of data from every individual. I actually use a script that obtains the cram file from here (ftp.1000genomes.ebi.ac.uk) and then processes it to obtain the fasta file.
The problem is that cram files are large, slow to download and slow to process. It takes days to get the samples.
Is there a better way to get these samples in fasta format?
The script already makes use of samtools to retrieve only the part of the bam file it needs but doesn't help much. Cram files are still gigabytes large for only a few megabytes of data that I need.
I have the same problem with data from the 1000 genomes project.

Could you please [edit] your question and clarify a bit? You say "I just need some megabytes of data from every individual", what individuals? Everyone in the 1000 genomes project? Then, you say "I have the same problem with data from the 1000 genomes project., what other data are you looking at? Also, what do you really need? Does it have to be from the 1000 genomes project? Do you just need X megabases of sequence from Y individuals and don't care where it comes from? Finally, do you want fasta or fastq?

You can download HGDP data in FASTQ format here: https://www.internationalgenome.org/data-portal/data-collection/hgdp


text:
Error in apply(dat[, batch == batch_level], 1, function(x) { : dim(X) must have a positive length
I am new here so forgive me if I don't provide all the details in the first go. There is a very similar question (how to solve “dim(X) must have a positive length” at running ComBat function in R) but the solution did not work for me.
I am getting Error in apply(dat[, batch == batch_level], 1, function(x) { :    dim(X) must have a positive length when I use the sva library's comBat() function.
It works when I use the following:
> dim(log_proteins_norm_sorted)
[1] 5415   66

> log_proteins_norm_sorted[1:5,1:5]
      1027     1075     1116     1119    11211
1 22.50513 20.64312 21.56182 22.01455 23.17497
2 21.67946 20.90994 21.61928 21.92025 21.98503
3 21.06812 20.52647 20.72076 20.79732 21.38957
4 18.95485 19.70027 18.39504 20.20218 19.79887
5 22.75966 22.32306 22.32936 22.39271 22.84204

> dim(samples_pheno)
[1] 66 24

> samples_pheno[1:3,1:3]
     Trial Number  Trial Sample
1027      222       A   1027
1075      223       A   1075
1116      224       A   1116

> proteins_combat = ComBat(dat= as.matrix(log_proteins_norm_sorted), batch=samples_pheno$Batch ,  par.prior = T)
Found10batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data

So as you can see, this works. However, when I subset some of the samples and reduce the sample number from 66 to 55 and run comBat, I get the following:
> dim(log_proteins_MRD_sorted)
[1] 5415   55
> log_proteins_MRD_sorted[1:5,1:5]
      1027     1075     1116     1119     1217
1 22.55492 20.69292 21.61162 22.06435 21.86900
2 21.72926 20.95974 21.66907 21.97005 21.35203
3 21.11791 20.57627 20.77055 20.84712 21.01896
4 19.00465 19.75006 18.44483 20.25198 18.95107
5 22.80946 22.37285 22.37915 22.44250 22.50524

> dim(samples_pheno3)
[1] 55 24
> samples_pheno3[1:3,1:3]
     Trial Number  Trial Sample
1027      222        A   1027
1075      223        A   1075
1116      224        A   1116

> proteins_combat_MRD = ComBat(dat= as.matrix(log_proteins_MRD_sorted), batch=samples_pheno3$Batch ,  par.prior = T)
Error in apply(dat[, batch == batch_level], 1, function(x) { :
  dim(X) must have a positive length

In case you wanted to know, I am working with the same original file for log_proteins_MRD_sorted and log_proteins_norm_sorted. I just removed the 11 samples from the 66 samples and then normalised the data again (hence slight difference in values between the two objects.
For samples_pheno3, I run this code:
> #Subset the phenotype data to those samples with only MRD data
> samples_pheno2 <- samples_pheno[-c(5,6,11,16,21, 26, 29,30,32, 45,56),]
> MRD_NA1 <- rownames(samples_pheno2)[!is.na(samples_pheno2$MRD)]
> samples_pheno3 <- samples_pheno2[MRD_NA1,]

Also this is out my data looks like for the batches:
> samples_pheno3$Batch
 [1] 7  5  6  9  7  7  7  7  8  8  8  8  9  9  9  9  10 10 10 10 11 11 11 2  7  2  2  2  9  3  3  3  3  3  3  3  7
[38] 10 3  4  10 4  4  8  4  4  4  5  5  5  8  11 5  5  5
Levels: 2 3 4 5 6 7 8 9 10 11

So just by subsetting the data to 55 samples, it doesn't seem to work for some reason. Is it the way I am subsetting the data which is confusing the code?
I would be grateful for your advice.

With the help of my very capable local bioinformatician, I have found a solution. Essentially, the problem was that when I reduced the sample size from 66 to 55, there was a batch number which only had 1 sample in it. When I removed that sample from the analysis, and ensured that each batch had at least 2 samples in it, the ComBat() then worked nicely!
Hope this helps future googlers.


text:
how to change the UMAP use in the dimplot and feature plot
I would like to know how to change the UMAP used in Dimplot and FeaturePlot from Seurat: how we can get the x-axis and the y-axis like UMAP-1 and UMAP-2 if I want to use UMAP-4 and UMAP-5. How to modify the code? It is easy to change the PC by using DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "PCA") but if I changed to reduction = "UMAP", I got the error
"Error in Embeddings(object = object[[reduction]])[cells, dims] : subscript out of bounds Calls: DimPlot Execution halted".

Does anyone have any idea how to solve this problem?

that is because UMAP-4 and UMAP-5 are not present in your object...head(object[['umap']]@cell.embeddings), or str(object[['umap']]@cell.embeddings)

Yes, the UMAP-4 and UMAP-5 are not present in my project. So what should I do now?

What is the rationale for using UMAP-4 and UMAP-5? UMAP is different from PCA because you specify a dimension (usually) and it tries to put all information in those two dimensions.

Try this:
library(Seurat)
pbmc_small <- RunUMAP(pbmc_small, dims = 1:10, reduction = "pca", n.components = 5)

names(pbmc_small[["umap"]])
#> [1] "UMAP_1" "UMAP_2" "UMAP_3" "UMAP_4" "UMAP_5"

# plot
DimPlot(pbmc_small, reduction = "umap", dims = c(4, 5))

This should work with Seurat version 3.1.5
packageVersion("Seurat")
#> [1] '3.1.5'

Change n.components / max.dims (depending on your seurat version) to 5 when you call runUMAP and you should have dim 4 and 5 available.

Hi Pallie, i have tried to use "SAN <- RunUMAP(SAN, dims = 1:10, max.dim = 5L, dims.keep = 5L)"  and there still are no UMAP-4 and UMAP-5.


text:
Parallelizing WGCNA k-means clustering and merging smaller clusters
WGCNA carefully takes advantage of server parallelization in several functions, and I can switch this on with:
enableWGCNAThreads(nThreads=28)

That speeds up some functions significantly, but the two bottleneck functions in the pipeline (in my experience) are part of the "Projective k-means" operation, primarily "..k-means clustering.." and "..merging smaller clusters..".
Unfortunately, parallelization does not work for those two functions.
Am I missing something?  How can I parallelize these functions?  Otherwise I'm wasting HPC resources and cluster time.
For reference, I'm working on a 28-processor HPC node with 96GB RAM, and I'm trying to run analysis on a 15,000-square correlation matrix.

The bad news is that, indeed, projectiveKMeans is not parallelized and I am not sure how much of it is (easily) parallelizable. The good news is that with 15k features (genes) and 96GB RAM you don't need to run preclustering at all. Just analyze the whole data set in one block.
If you for some reason insist on splitting the data into blocks, make sure the BLAS R uses on the cluster is a fast BLAS, ideally one that can parallelize matrix operations. A lot of time in merging of small clusters is spent calculating singular value decompositions, which becomes much faster using a good BLAS. HPC clusters will often have Intel KML available, or the admins could install OpenBLAS.

Thanks for each suggestion. Both helpful!  I definitely don't need to run preclustering, makes sense.  Will also make sure we are running Intel KML or OpenBLAS, I suspect this will already have been considered and optimized.

Will WGCNA k-means run using GPU? I have a hunch it might do.

Okay I checked and the answer is no. You're stuck, my bad. Reference https://bmcsystbiol.biomedcentral.com/articles/10.1186/s12918-017-0420-6


text:
What is a read count?
A very basic question: What exactly is a read count in shotgun sequencing and how is it calculated? I'm particularly interested in understanding the definition of the read count in the context of shotgun metagenomic sequencing though I assume the same concept applies to other types of shotgun sequencing as well. Surprisingly, a Google search did not results in any answers!
In the context of metagenomics, the taxonomic profiling of metagenomic data provides a list of microbial taxa and their "counts". These counts are often normalized to account for genome size or are normalized to convert counts to proportions. Understanding the concept of a read count may help me better understand all these.
Thanks!

To add to the answer from Bioathlete, "reads" in the traditional genomic sense are the string of nucleotide data inferred from the sequencer from a given fragment of genetic material (could be Illumina, Nanopore, PacBio, 10X, etc.).
So read count first applies for all of the "reads" defined above generated by the sequencer. From a NovaSeq, the total number of bases from a single flowcell is up to 3 trillion bases at 150bp each read for the highest throughput flowcell (so a read count of 3,000,000,000,000/150 = 20B total read count).
The read count for a single Nanopore flowcell is roughly (30 billion bases/avg read length), but the read lengths can be SUPER long (100kb compared to 150bp is substantial). Long read technologies are generally used for different analyses than short reads.
So all of the above reads would be defined as "raw reads", meaning they are reads not yet analysed in any way to be used for data analysis. So using illumina as an example, for a single Whole Genome sequencing sample we might have 800,000,000 raw reads. We then use an aligner like bwa mem to map the reads to a single reference genome (or in your case you would use a different aligner to map to multiple genomes). From here we still have the concept of raw reads, but we also now have the total number of mapped reads. A typical good illumina library for a single reference genome will have >90% mapping rate for a human sample.
We can then say we have 800,000,000 raw reads, and 720,000,000 mapped reads with a mapping to a reference genome of 90%. So we can compare "raw read count" to "mapped read count" to get an idea of the quality of our library.
In the case of metagenomics, you would consider your "raw read count" to still be the total number of reads you got off the sequencer, but the "mapped read count" would be the sum of the reads mapped to each individual reference genome to get a total quantity of "mapped reads" which could give you an indication of library quality. I'm not sure what a good metagenomics library "mapped read count" is, but I'm sure you can find that information.
Lastly, we want to get an idea of the depth of coverage for each mapped read to the genome. This will tell us how well our mapped reads actually covered the reference genome or genomes we are interested in.
I hope this was helpful!

Typically read count is the total number of reads going into the analysis. It could be based off single or multiple sequencing libraries.
Also it can be used to describe the number of reads that align to a region of the reference.  Depth or coverage are also terms used in this case.
You will need to be more specific in your use case for an accurate definition.

I modified my question to provide my specific details. I'm not sure if I understand what you mean by "the total number of reads going into the analysis".

"In the context of metagenomics, the taxonomic profiling of metagenomic data provides a list of microbial taxa and their "counts"" .  In this case it is the number of reads assigned to the taxa from the overall input sequencing data.


text:
What is the significance of enrichment factor regarding machine learning methods?
I'm fairly new to the field of bioinformatics and ran into a question while reading a paper I found on bioRxiv. The overall setting of the paper is using multi-task deep neural networks for kinase classification. The paper that I'm referring to is Kinome-Wide Activity Classification of Small Molecules by Deep Learning (Allen et al., 2019) for anyone who's curious.
I'll get straight into the specific figure that I'm having trouble understanding.

I know that questions like these are a bit discouraged on SE, but I'm having trouble on how to start regarding the interpretation of this figure.
Based on my understanding, the point seems to be that multitaks DNN's have the highest overall EF/EFmax ratio amounts, but what does this ratio signify?
I did some research and according to Wikipedia an "enrichment factor" in terms of bioinformatics is a concept "to measure the added value of a search tool over another one or over the homogeneous distribution in the genome population." Does this mean that multitask DNN's are "better" than the other models, due to its higher number of EF/EFmax points?
But if we look at the curves on the right, it seems to say that for fewer number of compounds, the ratio for multitask DNN's is lower...

I haven't read the paper, I don't do small molecules at all. The thing that is interesting is the multi-task DNN (deep neural network). Any chance you can explain its architecture please? DNN will possibly outperform the machine learn ing methods architecture and data set size depending. That is the result of this investigation, but more important is how they constructed the DNN, if they have borrowed it from natural language processing - I am not impressed. If they constructed it de novo that is a state-of-the-art approach.

The final question, given this is very likely under review are you reviewing the manuscript? The ratio is their accuracy estimate and is a standard approach.

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.

Ok I'll be brief, given the ms is likely under review.
Background in silico drug discovery has been part of deep learning for a long time 2015(?), and tensorflow was deployed. They worked out how to vectorise a small molecule and thereafter the entire field was.wide open to advanced ML.
Ratio The ratio is an accuracy index of the machine learning algorithms. The denominator is the total number of biological observations, the numerator is the predicted number. This 90% would be an accuracy index, i.e. the algorithm is correct 9 out of 10 times.
Result The plot shows, as the data set size increases the accuracy of the prediction increases most notably for the DNN and that conforms to all the theory. DNN/ANN only become powerful for 'big data', standard ML, i.e. random forest, naive Bayes, LR work well on small data sets, are fast and produce a classification which isn't wrapped up in a 3-4 Gig dump, thereafter performance saturates for routine ML methods.
The multi task DNN (MT-DNN) means that more than one loss function is being used, which is unusual in generic deep learning, but has been applied to natural language and drug discovery previously. Why this is being deployed is a long explanation, so I'll leave it there.


text:
Replicate papers (beginner)
I have a strong mathematical and computer background (also with machine learning) and would like to switch to bioinformatics. I know a fair amount of biology. But I find it hard to find real world practical training in this field. So I would like to try to replicate some bioinformatics papers (with their provided data). Good beginner papers for me would be shorter ones where the biological problem is well stated and the problem is also fitted in to the bigger scope (or is just easily understandable). Do you have any papers to recommend?

I'm sorry, but this really isn't a good place to ask this sort of broad, open ended question. We focus on specific problems here with specific answers. Asking for a list of papers isn't really on topic. You might have better luck on http://biostars.org

Bioinformatics is a big and very heterogeneous field with a lot of variation, so it's hard to recommend something that covers everything. That said, I think the best you'll do is to build familiarity with some of the common tools with some canned analyses.
One possible option would be some of the scRNA-seq workflows from the Pachter group, they seem to be fairly conscientious about documentation. Here is one such canned analysis, here is the paper, here is a blog post about it that suggests next steps and discusses at a high level.
scRNA-seq is just what's been trendy for a few years, next year it will be something different. But it at least gets you familiar with some of the file formats and concepts such as alignment, etc.
Hope that helps.

Thanks for your answer!
I'm not looking for something that covers bioinformatics, I'm looking for one or two bioinformatics papers which are replicate-able for someone that has started scratching the field. Most bioinformatics books mostly covers the different file-types and transformation between them, I want something with a little more substance but not to many layers (deep) in the biology part. e.g using machine learning to figure out the genes that are responsible for a particular genetic disease.

Well, then the linked analysis should help. It gets you from the raw data to data matrices, and points to clustering methods on those data, without getting too worried about the biology. Like many fields, that's what bioinformatics is, is composing sets of tools to get to datasets that you can plug into ML or whatever. There are many different ways of doing it, and many problems that can be studied by doing so. The linked analysis is one set of methods applied to one problem. Most of the "substance" happens upstream of the clean dataset, as with all analyses (even ignoring mechanism).


text:
Use of Electronic Phenotype in EHR
May I know what's the use of Electronic Phenotyping using EHR data?
I did refer this link but have few questions
I understand that Phenotype is a set of criteria that you apply on your EHR data to select patients of interest from EHR database.
However what I am trying to understand is the difference between ordinary SQL query and electronic Phenotype
If both are same,
1) what I would like to know is the use of electronic phenotyping? How it is used?
2) Why do we have pheKB repository and other repositories? What's the use in validating them in multiple sites?
3) Is it like identifying patients real time? I am not able to get the usefulness of this. I understand that we select patients of our interest (cohort), but i can just write an ordinary SQL and send it to my EHR database to identify patients. but what's special about Electronic Phenotype?
4) I also see that they use multiple conditions (Lab tests, Medications, Diagnosis etc) to identify patients. What am I missing here?
Can someone help me by giving an simple example? I am new to healthcare and trying to learn. Would really be helpful

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

No, a phenotype is a way of behavior or other observed characteristic of a person, resulting from a combination on its underlying biology and the environment it's in. Think things like, height, running speed, strength, etc.
In the context in which you're interested, this would be more properly, "the clinical presentation of a disease or disorder".

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

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

thanks for the response. Upvoted. Can you please help me understand this The point of an algorithm is to take the EHR data and progress to a diagnosis. ?

What is unclear about that?

I would suggest looking at the documentation for phenotips software and playing with the online demo.  This would show you some examples of phenotype data and how a database can be implemented to store it.  You may also want to familiarise yourself with HPO terms (Human Phenotype Ontology), a standardised ontology for recording human phenotypes (physical symptoms and characteristics which may be related to the underlying genotype).


text:
Finding differentially expressed genes between two samples

Hi, I was wondering whether we can find the differentially expressed genes between the Double-KO and the Shox2-KO. For the FindMarker() function, only see the differential expressed gene between different identities. I don't know whether we can find the differential expressed gene between only these two samples.

If I understand your question correctly, you can simply use SetIdent() to change the "default identity" to your samples and then use FindMarker() with the ident.1 = "Double-KO" and ident.2 = "Shox2-KO".
Alternatively, you can "collapse" or "summarize" your single cell data into a pseudo-bulk and then use one of the bulk RNA-seq packages.

Hi haci, that is a good idea. But I have tried to look the script of SetIdent(), find some difficulties to do it. Can you please give me more detail?

You can set the ident.use parameter of SetIdent() function to one of the columns in your meta.data slot of your Seurat object.

I'm not sure this is right.  SetIdent is used to assign new ident information; you give it a vector of cells to tag.  All you want to do is change the active ident that FindMarkers will use, and you do that with Ident().


text:
How could I match atom orders between a .mol2 and a .pdb?
When I was drawing conformation of Autodock result, I noticed that Autodock generated pdbqt from dlg has some structural problem (e.g. Benzene ring missing), and wants to correct this. However I noticed that the atom order of pdbqt doesn't match with the original mol2's, and caused some problems when I tried to fix.
So how could I reorder pdbqt according to mol2's atom coordinate?

Atom reordering is a common problem in compchemistry.
Rdkit (a python package) can do this, but it is limited by the formats it can read and mol2 files are a bit hit or miss. It works really well with SMILES, SMARTS and mol (sdf) files. But the writing may cause problems with Brookhaven pdb and mol2 files.
So the formats you have are both problematic. One is a pdbqt (a PDB file with partial charges (Q, generally Gasteiger calculated) and atom type (T) and the other a mol2. Let's pretend first that they are not.
Also the hydrogens might be made implicit depending on your mol2 file.
from rdkit import Chem
ref = Chem.MolFromMol2File('ligand.mol2')
target = Chem.MolFromPDBFile('ligand.pdb')
atoms = target.GetSubstructMatches(ref) #: List[int]
new = Chem.RenumberAtoms(target, atoms)
Chem.MolToMolFile(new, 'neoligand.mol')

Extra steps would be required for the protons and if your starting ligand was not from a file, but a SMILES.
mol = Chem.AddHs(mol) #protonate explicitly
Chem.GetSSSR(mol) #resonance fixing
AllChem.EmbedMolecule(mol) #initialise for 3d.
AllChem.UFFOptimizeMolecule(mol, maxIters=2000)

So this is all nice, but the formats may not work. Luckily openbabel can deal with this problem.
obabel -i mol2 input.mol2 -o pdbqt -O output.pdbqt

However, one thing to keep an eye out is what will be recalculated. Say, sdf files do not have partial charges or what PDB calls "atom names" and Rdkit calls "atom labels" (About the latter, I actually wrote a script to relabel atom labels in pdbs from mol files), while mol2 and PDBQT have.


text:
How to interpret the GROMACS MD Simulation RMSD Graph?
Hello brilliant researchers!
After virtual screening, I performed MD Simulation via GROMACS to check the stability of the protein-top hit complex by following this tutorial http://www.mdtutorials.com/gmx/complex/09_analysis.html
After that, I facing some difficulties in interpreting the resulted RMSD graph as it is quite different from most of the published RMSD graph.
Based on the graph attached, how to interpret it? Is AChE-dexketoprofen complex (Red) more stable (or conserved) than BuChE-dexketoprofen complex (Black)?

The opposite-ish.
Your graph is showing only the ligand, but RMSD is generally used for structures as a whole —so I'll mention about whole protein for future readers.
RMSD is a measure of how much does a structure deviate in position from a reference. Your ligand in the protein in red "wiggles" more than the one in black.
Structure RSMD
If this were RMSD of the protein (and not the index group for the ligand), it would not really mean that the black protein is more stable. Stability is to do with force-field energetics not position. Some protein have two state positions, say MBP is a flexing dumbbell. This would be seen as a jump in the RSMD graph, similar to that seen for the black trace.
However, the command you used was based on gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg from the page you linked. Key here is the index -n index.ndx, which limits it to the ligand. Before this happened gmx trjconv aligned the protein in the trajectory, so that common mistake was not done.
Ligand RSMD
In terms of ligand RMSD, lower deviation/less wobble correlates with more stable binding, but energetics (gmx energy) will tell you how strongly does it bind.
The black trace has a jump, so your ligand in the black trace (which wobbles less) has an alternative binding mode. This is why you are looking at RSMD.
To explore this further you can:

look at the trajectory
in PyMOL using the command split_states twice, once for the first and then the last state, followed by align of the two
If it were the whole protein, looking at the per bond matrix -bm option in Gromacs rms and spotting the change.

3Å??
I strongly suggest checking the trajectory because both have a deviation of 3Å, which for a ligand is not good.

Did the energy minimisation converge properly?
If it is because the proteins themselves wobble substantially, revisit gmx trjconv with a focus on the active site minus the ligand.
Is the ligand binding but leaving one substituent group to wobble?

Energetics
As with docking, the "score" is free energy difference (in kcal/mol). This is not the same as k_off and k_on and their ratio K_d (Dissociation constant). For those you need to do steered MD or have a way longer simulation (not recommended). The reason being is that the ligand has to make its way into the active site. However, you have a secondary binding mode in black, so have two free energy scores, this means that there is a stronger k_on.


text:
What is 'k' in sequencing?
When a DNA sequence is sequenced, I've only ever dealt with A,T,C,G and N which indicates un-identifiable bases. However, I came across a 'k' recently and I had asked another researcher who gave me an answer for what 'k' represents but I don't quite recall. It can't just be an anomaly in that one sequence file. Any ideas?

K is also the lysine amino acid (and k-mers are something else of course)

See IUPAC codes:

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

It is recommended you learn the degenerate nucleotide code. In sequencing it can signify poor quality sequence data, but in primer design it is useful. R (mutation within purines) and Y (mutation within pyrimidines) are common. K, a purine to pyrimadine or pyrimadine to purine mutation is, in my opinion, rare. I would treat a K mutation with caution and consider the triplet codon around it, i.e. if it is part of a protein gene.
Most phylogeneitcs programs will work with the degenerate nucleotide code, so in theory you can still obtain useful information with it.


text:
Most scalable pseudotime ordering algorithm
What algorithms for linear Pseudotime trajectory construction (diffusion-based) are the most scalable to large datasets?
I'm currently using Slingshot based on the recommendation in this manuscript: https://www.biorxiv.org/content/biorxiv/early/2018/03/05/276907.full.pdf
Slingshot is nice and easy to use, but it turns out it is terrible at scaling to big datasets.
For 2000 highest variable genes and a random sampling of a variable number of cells in my ~80% sparse expression matrix I get the following computation times in R (96 GB RAM):
#cells / Time (min)
100 / 0.7
1000 / 4.3
2000 / 12.3
5000 / 71.1

At this rate it would take nearly 1.5 days to process 200,000 cells and 4.5 months to process 2 million cells based on 2,000 HVG.
I feel like Slingshot is limited by a fundamentally unnecessary comparison of diffusion for all cells against all cells, rather than preclustering and then finer resolution of pseudotime across cluster edges.
What algorithms are better, but still implement a similar realization of pseudotime ordering?  I just need a single linear trajectory.

Well, in the absence of any comments, and being unable to find any scalable solution, I am now writing a function to granularly order cells at low resolution, then incrementally resolve high-resolution ordering. All using slingshot. Will post up when done.

Thank you! That would be appreciated!

As it turns out, pseudotime algorithms are exceedingly difficult to parallelize. Pseudotime does not scale linearly to a subset of data compared to the full dataset, and as such a course-grained approach to parallelizing pseudotime calculations on random bins of the data, followed by bin merging, does not work. These values will never converge on the actual full-dataset solution.
I put in a pretty good effort trying to parallelize slingshot, but with no success. I realize my code may be of limited utility without more explanation, but I'm not putting more time into this. FWIW:
This code:

Finds computationally optimal bin size for slingshot based on a number of genes and features
Randomly splits the SingleCellExperiment object into that number of bins and runs slingshot on each bin
Orients each bin so that a cell with a pseudotime of 1 in the first bin is similar to a pseudotime of 1 in every bin. Pseudotime can get flipped depending on the sampled subset.
Merges all bins into a single matrix so each cell has a assigned pseudotime based on randomly sampled bins.
Repeats steps #2 to #5 with a different random seed for binning every time.
Measures convergence towards actual solution at three pseudotime values, where convergence is defined as correlation of cell indices along the trajectory.

I abandoned the project after observing that more iterations to not change the convergence towards the actual solution. While pseudotime values close to 0 recapitulated the actual solution quite well, pseudotime values towards 100 failed to correlate at all with the actual solution. There was no convergence with more iterations.
Maybe this will be helpful to somebody, which is why I'm posting.
# Parallel Iterative Sampling with Slingshot
# pissshot is a course-grained estimator of the actual slingshot solution, and iteratively converges towards the actual solution
# pissshot randomly partitions cells in a large database into bins, bin size is determined by finding the computationally optimal number of cells
# pissshot is run independently on each bin, pseudotime vectors are aligned across all bins, and then all bins are merged
# pissshot repeats this process and the average of pseudotime assignments for each cell is taken at the end of each iteration.
# The convergence of the blingshot model at the end of each iteration is calculated by running slingshot on a small sample of adjacent cells at a random pseudotime interval and measuring correlation with the predicted model
# When satisfactory convergence is reached, blingshot returns a singlecellobject, just as slingshot would

pissshot <- function(sce, num_iterations = 20){
    cat("Running blingshot on", dim(sce)[2],"cells and",dim(sce)[1],"features\n")
    cat("  Step 1/5: Finding computationally optimal increment sizes\n   ")
#    increment_size <- FindOptimalIncrementSize(sce)
    increment_size <- 688
    num_increments <- ceiling(dim(sce)[2]/increment_size)-1
    cat("\n  ...Increment size of",increment_size,"cells is computationally optimal\n")
    cat("  Step 2/5: Course-grained pseudotime assignment across",num_increments,"increments")
    ps <- list()
    convergenceArr <- c()
    for(iteration in seq(from=1,to=num_iterations,by=1)){
        cat("\n   ... iteration",iteration,"...\n")
        cat("         randomly assigning cells to",num_increments,"bins\n         0%.")
        bins <- PartitionBins(sce, increment_size = increment_size, seed = 7*iteration)
        cat("100%\n         running slingshot on each bin\n         0%")
        cell_pseudotimes <- c()
        for(b in seq(from=1,to=length(bins),by=1)){
            bin_with_pseudotime <- suppressMessages(slingshot(bins[[b]]))
            # Find right cell and left cell in pseudotime array (cell with highest/lowest pseudotime values)
            pcells <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1)
            rownames(pcells) <- colnames(bin_with_pseudotime)
            pcells <- data.frame(lapply(pcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(pcells))
            pcells.sorted <- pcells[order(-pcells[,1]), , drop = FALSE]
            right_cell_ID <- rownames(pcells.sorted)[dim(pcells.sorted)[1]]
            left_cell_ID <- rownames(pcells.sorted)[1]
            right_cell_pos <- match(right_cell_ID,colnames(bins[[b]]))
            left_cell_pos <- match(left_cell_ID,colnames(bins[[b]]))
            right_cell_logcounts <- assays(bins[[b]])$logcounts[,right_cell_pos]
            left_cell_logcounts <- assays(bins[[b]])$logcounts[,left_cell_pos]
            if(b==1&&iteration==1){
                ref_right_cell_logcounts <- right_cell_logcounts
                ref_left_cell_logcounts <- left_cell_logcounts
            }
            # Figure out whether this pseudotime in this bin has flipped relative to reference (first bin), and if so, flip pseudotime values
            cor_opposite_ends <- mean(cor(right_cell_logcounts,ref_left_cell_logcounts),cor(left_cell_logcounts,ref_right_cell_logcounts))
            cor_same_ends <- mean(cor(right_cell_logcounts,ref_right_cell_logcounts),cor(left_cell_logcounts,ref_left_cell_logcounts))
            flipped <- FALSE
            if(cor_opposite_ends > cor_same_ends){
                flipped <- TRUE
            }
            if(flipped==TRUE){
                # find max pseudotime in colData(bin_with_pseudotime)$slingPseudotime
                # Recalculate pseudotime in bin_with_pseudotime as max-value
                pseudotimes <- colData(bin_with_pseudotime)$slingPseudotime_1
                colData(bin_with_pseudotime)$slingPseudotime_1 <- max(pseudotimes)-pseudotimes
            }
            cells_with_pseudotime <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1,colnames(bin_with_pseudotime))
            cell_pseudotimes <- rbind(cell_pseudotimes, cells_with_pseudotime)
            cat(".")
        }
        cat("100%\n")

        ps[[iteration]] <- cell_pseudotimes[order(row.names(cell_pseudotimes)), , drop = FALSE]

        pss <- data.frame(ps)
        if(iteration>1){
            convergence10 <- measureConvergence(pss, sce = sce, ptime = 10, increment_size = increment_size)
            convergence50 <- measureConvergence(pss, sce = sce, ptime = 50, increment_size = increment_size)
            convergence99 <- measureConvergence(pss, sce = sce, ptime = 99, increment_size = increment_size)
            cat("   convergence at 10:",convergence10,"... at 50:",convergence50,"...at 100:",convergence99,"\n")
            convergenceArr <- rbind(convergenceArr,c(iteration,convergence10,convergence50,convergence99))
        }
        # Iterate until the model converges towards the actual slingshot solution, within an indicated fraction (i.e. convergence = 0.02)
        # The convergence measure compares how similar the model is to the actual slingshot solution, i.e. a convergence of 0.02 means that cell pseudotime values are within 2% of the actual solution
        # Assess this by selecting a pseudotime at random, analyzing the following increment_size cells, and comparing accuracy of pseudotime measurement across that trajectory with the composite of the course-grained analysis

        # need to study how convergence changes over pseudotime
    }

    # remove cells with really high standard deviations

    return(pss)
#    return(convergenceArr)
}

measureConvergence <- function(mat, sce = sce, ptime = 45, increment_size = 500){
    # this line is throwing an error
    mat[,"avg"] <- apply(mat[,1:dim(mat)[2]],1,mean)
    mat <- mat[order(mat[,"avg"]), , drop = FALSE]
    # select a random pseudotime increment
    mat <- mat[mat[, "avg"] > ptime,]
    mat <- mat[1:increment_size,]
    cell_IDs <- rownames(mat)
    # pull out this list of cell_IDs from the sce object and run slingshot
    sce2 <- suppressMessages(slingshot(sce[,cell_IDs]))
    cpcells <- as.data.frame(colData(sce2)$slingPseudotime_1)
    rownames(cpcells) <- colnames(sce2)
    cpcells <- data.frame(lapply(cpcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(cpcells))
    # sort cpcells by increasing pseudotime
    cpcells <- cpcells[order(cpcells[,1]), , drop = FALSE]
    ordered_cells <- cbind(rownames(mat),rownames(cpcells))
    # get the index of each cell in the dataframe, run a correlation on how well the indices line up
    cell_indices <- c()
    for(cell in cell_IDs){
        new_indices <- c(which(ordered_cells[,1] == cell), which(ordered_cells[,2] == cell))
        cell_indices <- rbind(cell_indices,new_indices)
    }
    convergence_measure <- cor(cell_indices[,1],cell_indices[,2])
    plot(cell_indices[,1],cell_indices[,2])
    return(convergence_measure)
}

PartitionBins <- function(sce, increment_size = 500, seed = 123){
    bins <- c()
    while(dim(sce)[2]>=increment_size){
        if(dim(sce)[2]<increment_size*2){
            bins <- c(bins, sce)
            sce <- sce[,1]
        }
        else{
            set.seed(seed)
            train_ind <- sample(seq_len(ncol(sce)), size = increment_size)
            bins <- c(bins, sce[,train_ind])
            sce <- sce[,-train_ind]
            cat(".")
        }
    }
    return(bins)
}

FindOptimalIncrementSize <- function(sce){
    # Find the computationally optimal increment_size for grainy alignments using 5 increment points (100, 500, 1000, 2500)
    res <- c(100,250,500,750,1000,1500,2500,5000)
    runtimes <- c()
    for(i in res){
        runtimes <- c(runtimes, suppressMessages(system.time(slingshot(sce[,0:i])))["elapsed"])
    }
    mat <- as.data.frame(t(rbind(res,runtimes)))
    c <- summary(lm(mat$runtimes ~ poly(mat$res, 2, raw=TRUE)))$coefficients[,"Estimate"]
    predicted_runtimes <- c()
    for(i in seq(from=100, to=5000, by=1)){
        predicted_runtime <- (c[1] + c[2]*i + c[3]*i^2)*(dim(sce)[2]/i)
        predicted_runtimes <- rbind(predicted_runtimes,c(i,predicted_runtime))
        cat(".")
    }
    increment_size <- predicted_runtimes[which.min(predicted_runtimes[,2]),1]
    return(increment_size)
}

There is a recent publication that compares methods for trajectory inference. In Figure 2 in terms of scalability PAGA, PAGA Tree and Component 1 get the highest rating. The authors offer a shiny app to help select the best method for your dataset.


text:
How to Convert From Gene Location to SNP Location
I am new to this. I have an SNP array composed of digits [0,1,2] for chromosome 21. How do I find the start:end locations for the region that is "21q22.11" step by step?

It seems your title doesn't really match to the content of your question, could you fix that?

You first need to get the coordinates of the "cytoband" 21q22.11, which I would search for using the UCSC table browser. You need the following settings:
Clade: Mammal
genome: Human
assembly:
group: Mapping and Sequencing
track: Chromosome Band
table: cytoBand
region: position chr21
output format: BED
(you can leave other options as default)
That will give you a bed file of cytobands on chr21. You'll need to search (e.g. grep) for your cytoband: q22.11.


text:
Downloading dataset from SRA (SOLiD Platform)
Trying to download colorspace data from SRA, but getting an error
abi-dump -A SRR1657115.sra
abi-dump.2.9.6 err: item not found while constructing within virtual database module - the path 'SRR1657115.sra' cannot be opened as database or table.

You have to use the accession name without the .sra suffix.
abi-dump -A SRR1657115


text:
Wrong format .ped in PLINK
I know that I don't have the correct format in my .ped file for PLINK but I don't know how to convert it.
My actual file is something like:
BSCAN BSCAN20016 0 0 2 0 0 1 1 1....
BSCAN BSCAN20017 0 0 2 0 0 1 1 1....
BSCAN BSCAN20018 0 0 2 0 0 1 1 1....

Where column 1 is FamilyID, column 2 is ID, column 3 is FatherID, column 4 MotherID, column 5 Sex and then the genotype of each individual.
I only need to create binary files, so I tried to run:
plink-1.9-rc --file admixed --make-bed --no-pheno --cow --out admixed

but the program shows: Error: half-missing call in .ped file at variant 1, line 1.
Anyone know how can I change the format or recode it to make a correct input?
Thanks in advance!

Please consider the response below - its a good response - and consider upvoting/accepting.

So, columns 7 and 8 should be the genotype of variant 1 (each column is one of the alleles). You have them as 0 1, which means absent and major allele.
I think you're just missing a column. From the doc:
Column 6 is Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)
Ideally you would know which SNPs are heterozygous based on your SNP calling and have such SNPs as 1 2, but if you have a homozygous organism or if you're working with data that's been flattened you can simply generate the .ped file with the same number twice for each SNP.

Thanks for your reply. So, I need to phase the data? Because now I have one column per SNP.  I don't have the phenotype column but I specify in the command "--no-pheno"

Yeah you're right that --no-pheno should mean it can read them fine, but it says specifically that the first variant in the first column is half-empty. Can you try after adding a phenotype column? You can use -9 to indicate you don't have one.

Oh, is your example not your real data maybe? I assumed it was real aside from there being more lines and columns

Actually, it is quite similar to my real data... but I have 5564 individuals (rows) and 46465 SNPs (columns). I will try adding the phenotype column...

I have the same error with the phenotype column.

On the first line, if you do a cat $file | head -n 1 | awk '{print $6,$7}' do you get 0 1? (or 1 0?). Because that's what I think the error message is saying you have

In my real data, next to the column of phenotypes, I have 0 0 1 0 1 1 1

I think that I found the main problem. I should phase the data.

No, .ped files do not contain phase information, but they do require two columns per variant, and those columns should contain allele codes, not alternate-allele counts.
So, if your first variant is an A/G SNP and your second variant is G/T, columns 6 and 7 of your file should always be "A A", "A G", "G G", or "0 0"; and columns 8 and 9 should always be "G G", "G T", "T T", or "0 0".


text:
Searching for proteins by primary and secondary structure
I need to find peptides/proteins with a certain sequence of amino acids in a certain type of secondary structure. For example, "ALA PHE GLY" in alpha-helix. Is there any tool that does that?

Your request is rather specific, so there is no tool...
It is a straightforwards request though. I'd impose the requirement that the structures be solved and no predicted secondary structure. You can do it two ways if using Python.
PDB
After downloading all PDB structures, in python with PyMOL (installable with conda, tricky without) you'd do something like
import pymol2, os
from multiprocessing import Pool

def find_motif(file):
    motif = 'AFG'
    try:
        with pymol2.PyMOL() as pymol
            pymol.cmd.load(file)
            seq = pymol.cmd.get_fastastr()
            if motif not in seq:
                return None
            d = {x: []}
            # str.find will cause issues if there is a gap also structure may not start at one.
            pymol.cmd.iterate('resn ALA+PHE+GLY', 'x.append((resi, resn, ss))', space=d)
            sequence = -1
            i = 0
            for residue in d['x']:
                if residue['resn'] = 'ALA':
                    sequence = 0 # lazy enum
                    i = int(residue['resi'])
                elif sequence != -1 :
                    if int(residue['resi']) == i+1 and residue['resn'] = 'PHE':
                        sequence = 1
                    elif int(residue['resi']) == i+2 and residue['resn'] = 'GLY':
                        return (file, i)
                    else:
                        sequence = -1
    finally:
        return (file, None)

if __name__ == '__main__':
    p = Pool(12)
    structure_files = [f for f in os.listdir('pdb') if os.path.splitext(f)[1] == '.pdb']
    verdicts = p.map(finder, structure_files)
    print([(f,v) for f,v in verdicts if v is not None])

The end result will be a list of PDB with where in the sequence (in that PDB code).
(I wrote this in the reply box, so there may be a typo or two).
PyMOL has a selection algebra called pepseq, which matches missing atoms, but I was a whimp and played it safe by being more verbose. Also, a for loop can be used instead of an iterate. But I think it is slower.
atoms = pymol.cmd.get_model(f"pepseq {target} and name CA")
for atom in atoms.atom:
      if atom.ss == ' ...

Uniprot
In Uniprot, if a protein is crystallised, there will be the feature entries corresponding to the secondary structure. It is sheet or helix. No turns, pi-helices or 3_10 helices. There is no Biopython method for parsing Uniprot XML. I wrote myself a parser and it is not too horrid thing to do —although I monkeypatched elementtree. Alternatively, xmlschema.XMLSchema.to_dict can be used to convert the XML to a dictionary. The advantage of Uniprot is that the residues ids will not require correcting. The disadvantage is parsing XML.
My parser is specific for my application and it is poorly documented (so would not recommend it) but but... I can share it and human data. I'd be used like:
 import os
 from michelanglo_protein import ProteinCore
 ProteinCore.settings.startup(data_folder='protein-data')
 p = ProteinCore().gload(file=os.path.join('gpickle','Q86V25.gpz'))
 #print(p.asdict())
 if 'helix' in p.features:
        n = p.sequence.find( ...
        for ... p.features['helix'] #('helix', 'turn', 'strand')
        if x < n < y ...

Thank you. Parsing it on one's own doesn't look complicated. Still, I somehow balk at the idea that my needs are so specific... You are positive there's no such tool?

Proving a negative is hard, so I can only say I have never read of one.But I had a similar problem some time back and did it manually. Furthermore, logistically a server to do this would need its data pre-generated but then the combinatorial issues will be huge so would not work...

there is msdmotif https://www.ebi.ac.uk/pdbe/docs/Tutorials/msdmotif/introduction.html but I suspect the query might be a bit slow as there would be a lot of hits

I just realised that in my comment I said something incredibly daft. A server is totally possible, they just need to store sequence+SS type sequence as opposed to all combinations. Apologies.


text:
Removing rows containing some strings
I have a big data frame
> dim(anno_maf)
[1] 449892    132
>

In chromosome column I have chr 1 to 22 and another contigs like
> unique(anno_maf$Chromosome)
 [1] "chr1"                  "chr6"
 [3] "chr16"                 "chr17"
 [5] "chr20"                 "chr2"
 [7] "chr3"                  "chr4"
 [9] "chr14"                 "chr19"
[11] "chr5"                  "chr10"
[13] "chr9"                  "chr12"
[15] "chr13"                 "chr11"
[17] "chr22"                 "chr7"
[19] "chr15"                 "chr8"
[21] "chr18"                 "chr21"
[23] "chrX"                  "chrY"
[25] "chr4_gl000194_random"  "chr17_gl000205_random"
[27] "chrUn_gl000241"        "hs37d5"
[29] "chrUn_gl000219"        "chrUn_gl000234"
[31] "chr1_gl000191_random"  "chrUn_gl000211"
[33] "chrUn_gl000224"        "chrUn_gl000225"
[35] "chr17_gl000203_random" "chrUn_gl000212"
[37] "chrUn_gl000243"        "chrUn_gl000214"
[39] "chrM"                  "chr1_gl000192_random"
[41] "chr7_gl000195_random"  "chrUn_gl000232"
[43] "chr4_gl000193_random"  "chr19_gl000208_random"
[45] "chrUn_gl000226"        "chrUn_gl000218"
[47] "chr9_gl000199_random"  "chrUn_gl000217"
[49] "chrUn_gl000229"        "chrUn_gl000216"
[51] "chrUn_gl000231"        "chr9_gl000198_random"
[53] "chr17_gl000204_random" "chrUn_gl000220"
[55] "chrUn_gl000235"        "chr11_gl000202_random"
[57] "chrUn_gl000222"        "chrUn_gl000240"
[59] "chrUn_gl000233"        "chrUn_gl000230"
[61] "chrUn_gl000213"        "chrUn_gl000238"
[63] "chr19_gl000209_random" "chrUn_gl000237"
[65] "#CHROM"

  How I can only keep rows having chr1 to chr22 and removing the rest?

You can create a vector consisting of values that you would like to keep (or discard):
to_keep <- paste0("chr", seq(1,22))

> to_keep
 [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"  "chr10" "chr11"
[12] "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18" "chr19" "chr20" "chr21" "chr22"

your_df[your_df$Chromosome %in% to_keep,]

The data.table way, which should be faster for large data, would be:
your_df[Chromosome %in% to_keep]

I don't know what your goal is but you might want to keep chromosomes X and Y as well.


text:
Merging two data frames
I have a data frame like with more samples like
> head(a[1:2,1:4])
             ID               ids pre.or.post.treatment
1 OCCAMS/AH/120 LP6005409-DNA_F01                   pre
2 OCCAMS/AH/126 LP2000325-DNA_A01                   pre
            TRG
1 non-responder
2 non-responder
> dim(a)
[1] 160   4
>

Another data frame with part of my samples but more characteristics for each sample like
> dim(b)
[1] 60 50

> head(b[1:2,1:6])
             ID       TRG pre.or.post.treatment          ID.1
1 OCCAMS/AH/212 responder                   pre OCCAMS/AH/212
2 OCCAMS/AH/279 responder                   pre OCCAMS/AH/279
  DiagnosisDate FE.LastSeenDate
1      8/8/2013        6/6/2014
2     6/10/2014       11/8/2018
>

How I can  merge these data frame getting a new data frame with my all 160 samples columns of both a and b

Assuming your column ID in both dataframes a and b shared same IDs, and as your dataframe are not particularly heavy, you can have the use of the function left_join from dplyr package as this:
library(dplyr)
DF <- left_join(a,b, by = "ID")

@dc37's dplyr solution is quite elegant, but you can also use base R's merge:
merge(x = df_a, y = df_b, by = "ID", all.x = TRUE, all.y = TRUE)

You can tweak various parameters, such as use by.x and by.y instead of just by, change all.x and all.y based on requirements to get your results.
This is a pretty common R question by the way, and if you had googled "merging data frames", you'd see a lot of options.


text:
how can i calculate p value of miRNA raw reads counts using R?
working on lung cancer project and need to find out which genes are deferentially expressed and which are not. for that i a have downloaded miRNA data set from GEO id is GSE110907. And the file i got has to variable gene_id and read counts. from that i have to derive p value from that data/file.

For context, this question was already asked at biostars https://www.biostars.org/p/417111/, with some follow-up comments.

I would start by re-formulating the question in this post: Your end-goal seems to the differential expression analysis (DEA) of your samples (and you have to have to levels for this like treatment-control, mutant-WT, etc ...). P-values are context-dependent and should be calculated based on your research question/hypothesis.
DEA has been addressed by many researchers and R has a lot of packages to help you with this; DESeq2, edgeR, limma, etc ...
I would go to the Biconductor website, search for one of the three above, and go through their rich vignettes and would come back here for specific questions.