Why does this human bam file only have one copy of each chromosome?

The maternal and paternal copies of a chromosome are called haplotypes. Many metazoans (animals) are diploid and have maternal and paternal chromosome contribution during sexual reproduction, not just humans as your question states.

Response to Q1

Your question, in other words, is: Why do .bam files not differentiate between haplotypes?

Your question gets to a more fundamental kernel of how most people do "genomics" today. When most people assemble a genome or create a reference genome, they actually just produce a fasta file that contains a single sequence for each chromosome. This, of course, is not biologically accurate as there are two very distinct sequences per chromosome. Most genomes only report a chimeric consensus between the two and call this a reference. These reference genomes are haploid, or haplotype-collapsed.

This is where your question comes in: You have mapped transcripts that biologically are derived from two different haplotypes to a reference genome that only contains one consensus sequence that (inadequately) represents both haplotypes. As a result, at a single locus in the bam file there will be mapped reads from both haplotypes. If your reference fasta file contained the exact sequence of both haplotypes, if you mapped reads to the reference, and if you only look at primary alignments, the reads will mostly map to their haplotype-of-origin.

This gets into another topic called phasing, wherein the order and orientation of most of the polymorphisms unique to each haplotype are determined using sequencing data. There are some problems with this as it relies on properly detecting variant sites. Software like GATK and others can find single nucleotide polymorphisms (SNPs) if there is good sequencing coverage, however detecting insertions and deletions is much more difficult. This gives a very SNP-biased view of the haplotype differences in any given genome. After finding variant sites, phasing software like  hapcut2 determines which variants fall on which haplotype and outputs blocks of variants that are thought to belong to the same haplotype.

Phasing alone is not enough to accurately reconstruct the exact sequence of both haplotypes due to the inability to detect all variants with read mapping. The gold standard of the future is diploid de novo genome assembly, in which each haplotype is assembled independently. This is an active area of research for people who develop genome assemblers and is tightly linked to the advances in PacBio and Oxford Nanopore long reads. This paper on trio canu is a good start for learning about one successful dioploid assembly technique.

Response to Q2

If you want to check ploidy of all of the chromosomes then you need at least 10X coverage of whole-genome shotgun data and a tool like smudgeplot and genomescope.

If you are trying to check local duplications or whole-chromosome duplications for something like a cancer sample you would also need whole-genome shotgun data. WES data does not give reliable information on ploidy since the read count at a given locus is dependent on the transcription level, not the actual quantity of homologous chromosomal DNA for the region in question.

If you are trying to (use a bam file to) look for evidence of chromosomal duplication or duplication of a locus in cancer samples you would look for an increase in whole genome shotgun coverage at that locus compared to a known normal locus. For example if all chromosomes but one have a coverage of 22, and one chromosome has an average coverage of 33, this is evidence of a trisomy. This logic can be applied to smaller regions as well (barring repetitive regions, paralogs, et cetera).


Why aren't clustal alignments stable under deletion of some of the sequences?

As the graph of sequences depends on the sequences in the input, if you change the input sequences (by deleting some sequences) you change the graph.
The alignment between the sequences changes, and thus the relative distance between them changes too.

The deletion of some sequences won't cause any change on other sequences only if the full space of the C(n 2) pairwise comparisons would be done. However for a large number (and long) of sequences, this would require too much computational effort and time, making such methods as the k-tuple distance competitive.


Is re-normalization of RNAseq data recommended for analysis of gene subsets?

FPKM, RPKM are ways of normalising for both sequencing depth and gene length. You first calculate a scaling factor (SF = #read/1Mil),you then divide the genes counts for the SF(RPM = # counts/SF), and finally you divide by the length of the gene (RPM/length(gene)). Each sample has a different SF meaning that you should not compare values from different samples.

To fix this, you have to normalise across samples, by upper quartile normalise using a subset of genes. In this case, the SF is calculated based on the upper quartile (70%) of this shared and immutable subset, and it is therefore the same for all the samples.

If you subset your genes, and you are using the normalised values for downstream analyses, then you need to re-scale your values.

I personally try to avoid FPKM,RPKM,TPM (if I can) and always start from raw counts.

I would suggest to have a look at the documentation of packages like DESeq2 and edgeR for detailed and better explanations.


Check if a string is ensembl gene id or not in python

Fixed it myself actually, not that hard. I found regular expression for ensembl ids here:

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

And used the re package in python:

import re

ensembl_regexp = 'ENS[A-Z]+[0-9]{11}'

is_ensembl = np.all([bool(re.match(ensembl_regexp, gene)) for gene in gene_symbols])
    if is_ensembl:
        gene_symbols = list(f['ensembl_ids'])

Where gene_symbols is a numpy array that I need to check whether it has ensembl ids, and if yes, switch to using gene symbols which are going to be under ensembl_ids key due to a mistake of someone who made the hdf5 dataset.


Why DoHeatmap Does not show all genes in genes.use?

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

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


Filtering Sequences (entries) by headers ID from a Fasta file database

Using Biopython:
from Bio import SeqIO

for record in SeqIO.parse("input.fa", "fasta"):
    if record.description.split()[-2:] == ["Pseudomonas", "aeruginosa"]:
        print(record.format("fasta"))


Biplot scores from canonical correspondence analysis

Canonical correlation analysis (CCA) not correspondence canonical analysis. Biplots are the same as in PCA: the representation in a new space of the samples and the relative position of the variables on the new space.

However, CCA seeks the values of one block that maximize the correlation with the other block. Here you don't have an independent variable that you maximize from the variables you introduced on the CCA. So the interpretation is much trickier as there isn't a direction: WaterC and FallenTwigs show the opposite effect than BarSand and CoverMoss, but there is a group of variables that explain what is the effect of water on the plants(?)


move files of a different group to separate directory

You need to bring together the information in the two tables, I prefer R for such actions. The common columns between the two tables do not quite contain the same information hence additional data cleaning was required. Moreover the table with treatment information had to be changed to long format with melt. And I have filtered the merged_data so as to conation only "pretreatment" info but you can add a if/else statement to the for loop to control how/where you will move your files.

library(data.table)

clin_name <- fread("clin_name.txt", header = TRUE)
clin_name$sample_id <- sub("/E/T.", "", clin_name$sample_id)

pre_post <- fread("pre_post.list", header = TRUE)
print(pre_post)
pre_post <- melt(pre_post,
     measure.vars = c("pre-treated", "post-treated"),
     variable.name = "treatment",
     value.name = "sample_id")

> print(pre_post)
      treatment sample_id
1:  pre-treated OC/AH/120
2:  pre-treated OC/AH/126
3:  pre-treated OC/AH/136
4:  pre-treated OC/AH/160
5: post-treated OC/AH/065
6: post-treated OC/AH/069
7: post-treated OC/AH/074
8: post-treated OC/AH/084

merged_data <- merge(clin_name, pre_post)
print(merged_data)

pre <- merged_data[merged_data$treatment == "pre-treated", "file_name"]

> print(merged_data)
   sample_id         file_name   treatment
1: OC/AH/120 LP6005409-DNA_F01 pre-treated
2: OC/AH/126 LP2000325-DNA_A01 pre-treated
3: OC/AH/136 LP6005409-DNA_E02 pre-treated

# invoking a terminal command to move each file to a directory called "pre"
for (i in pre$file_name) {
  command <- paste0("mv ", i, " pre/")
  system(command)
}


Automating sequence download from Wormbase

Further inspecting the source code of the downloaded page, I found the following:

    <span class="ui-icon ui-icon-close module-close" wname="sequences"

    Sequences

This suggested me to have a look at http://www.wormbase.org/rest/widget/transcript/K06C4.12/sequences and I obtained a page containing the elements found in the DOM inspector.

This can be parsed, for instance using the bs4 python module:

# To get the page
import requests
# To parse the page content
from bs4 import BeautifulSoup
# To get this into Biopython
from io import StringIO
from Bio import SeqIO

rest_template = "http://www.wormbase.org/rest/widget/transcript/{accession}/sequences"
rest = rest_template.format(accession="K06C4.12")
# This did not work with "html.parser"
soup = BeautifulSoup(requests.get(rest).text, "html5lib")

def generate_fasta(soup):
    name = None
    for div in soup.find_all("div"):
        if "ui-helper-hidden" in div["class"]:
            if div["id"] == "content":
                fasta = SeqIO.read(StringIO(div.text), "fasta")
                description_fields = fasta.description.split()
                if len(description_fields) == 4:
                    fasta.id = "_".join(description_fields[:2])
                    if name is None:
                        name = fasta.name
                elif len(description_fields) == 2:
                    # This assumes that the "conceptual_translation"
                    # is after the "unspliced + UTR" or "spliced + UTR"
                    # Otherwise, one would need to pass
                    # the accession as function parameter
                    fasta.id = f"{name}_translation"
                else:
                    fasta.id = "_".join(description_fields)
                yield (fasta.id, fasta)

records_dict = dict(generate_fasta(soup))
print(records_dict['K06C4.12_spliced'])
print(records_dict['K06C4.12_translation'])

Output:

ID: K06C4.12_spliced
Name: K06C4.12
Description: K06C4.12 spliced + UTR
Number of features: 0
Seq('ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTC...TAA', SingleLetterAlphabet())
ID: K06C4.12_translation
Name: conceptual
Description: conceptual translation
Number of features: 0
Seq('MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSS...SK*', SingleLetterAlphabet())


Are there any RepBase alternatives for genome-wide repeat element annotations?

Dfam has recently launched a sister resource, Dfam_consensus, whose stated aim is to replace RepBase. From the annoucement:

  Dfam_consensus provides an open framework for the community to store both seed alignments (multiple alignments of instances for a given family) and the corresponding consensus sequence model.

Both RepeatMasker and RepeatModeler have been updated to support Dfam_consensus.

I haven’t tried it yet but it looks promising.


Branch-models are currently not implemented in godon. In the PAML manual Ziheng Yang writes:

This model is very parameter-rich and its use is discouraged.

This was one of the reasons I did not implement it.
Disclaimer: I am the author of godon.


You can read in the normalized count table and don't normalize the data, but my advice here is not to do that. The rounding of the normalized matrix introduces some noise, but I think the larger issue is how sure are you that the table you are working with, is exactly a count table of normalized counts from DESeq2 ?
It really doesn't take much time to align and get count table again, rather than going back and forth wondering whether you did the normalization correctly.
I can show below with an example how it will look like, first we make an example dataset and obtain the normalized counts:
library(DESeq2)
set.seed(111)
sz = runif(6,min=0.5,max=1.5)
x = makeExampleDESeqDataSet(sizeFactors=sz,m=6)
x = estimateSizeFactors(x)
ncounts = counts(x,normalize=TRUE)
x = estimateDispersions(x)

Now suppose we feed this into a new DESeq2 object, but before that you need to do is check most of the size factors are close to 1:
estimateSizeFactorsForMatrix(round(ncounts))
 sample1  sample2  sample3  sample4  sample5  sample6
1.039179 1.045991 1.048276 1.045872 1.032420 1.052006

If you see something above that suggest a deviation from 1, something is wrong from the normalization and do not even think about using it.
Below you can use the normalized counts and you can see the difference is not huge:
dds = DESeqDataSetFromMatrix(round(ncounts),colData(x),~condition)
sizeFactors(dds) = 1
dds = estimateDispersions(dds)

Compare the new count table:
head(counts(dds,normalize=TRUE))
      sample1 sample2 sample3 sample4 sample5 sample6
gene1       0       2       0       1       0       0
gene2       9       8      16       6       1       2
gene3      23      40       8      38      11      41
gene4       0       0       0       6       8       0
gene5       3       2       0       5       1       2
gene6       0       1       3      19       4       9

And your input:
head(ncounts)
        sample1    sample2   sample3    sample4   sample5   sample6
gene1  0.000000  2.2658464  0.000000  0.9310803  0.000000  0.000000
gene2  8.693480  8.3081033 16.251453  5.5864817  1.121781  2.086619
gene3 23.472395 40.0299524  7.584012 38.1742919 11.217808 40.689076
gene4  0.000000  0.0000000  0.000000  5.5864817  7.852465  0.000000
gene5  3.477392  1.5105642  0.000000  4.6554014  1.121781  2.086619
gene6  0.000000  0.7552821  3.250291 18.6216058  4.487123  9.389787

We can check the dispersion, for lowly expressed genes you see some difference, but on the whole difference
with(elementMetadata(x),plot(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#FFC93C80"))
with(elementMetadata(dds),points(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#DBF6E980"))

Again, I would say I agree with @ATpoint comments that if you can, start from the fastq again.. don't work with data you have doubts about.


Find specific functional genes from a list?

BioMart. Filter by your list of genes and a GO term linked to acetyl transferase activity.


Interesting question I would recommend the effective population size $N_{e}$,  and error bars around this. Ne can be calculated via $F_{st}$ and is perhaps the most important parameter to determine the genetic population size.
Under the Fisher definition $N_{e}$ is the number of breeding individuals in the population.
The central issue with the approach you are suggesting is that population size does not equate to genetic population size. Thus the population could be huge, but if it is genetically homogeneous it will appear very small. In contrast the number of individuals could be very small, but the genetic variation between between them very large so $N_{e}$ would be huge.


Help with setting DimPlot UMAP output into a 2x3 grid in Seurat

how to get this to arrange the UMAP plots into a 2x3 grid

You can use the ncol parameter, which will set the number of columns in the "combined plot".
And for future reference, if you get the help page of the function you are using with ?DimPlot, you will see the .... This means that you can feed additional arguments to other functions that DimPlot() uses. Apparently,  CombinePlots() is used under the hood and has an ncol parameter.


Clustering information saved in Seurat object

I am not sure of the output of PrintFindClustersParams, as I haven't use it before, but try with object@commands. It returns a list containing every command and parameters used.
For example
...
$FindNeighbors.SCT.pca
Command: FindNeighbors(seurObj, dims = 1:20, verbose = vrb)
Time: 2019-12-11 12:41:12
reduction : pca
dims : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
assay : SCT
k.param : 20
compute.SNN : TRUE
prune.SNN : 0.06666667
nn.method : rann
annoy.metric : euclidean
nn.eps : 0
verbose : FALSE
force.recalc : FALSE
do.plot : FALSE
graph.name : SCT_nn SCT_snn

$FindClusters
Command: FindClusters(seurObj, verbose = vrb)
Time: 2019-12-11 12:41:17
graph.name : SCT_snn
modularity.fxn : 1
resolution : 0.8
algorithm : 1
n.start : 10
n.iter : 10
random.seed : 0
group.singletons : TRUE
verbose : FALSE
...

Re: update
I think it depends on how you built the object. FindNeighbors.SCT.pca should be result of running sctransform normalization. If you use the standard way, you may have other names (e.g.FindNeighbors.pca).


Is there a way to customize what observations are plotted in ComplexHeatmap?

What is fed into ComplexHeatmap is what it outputs! So if one of your genes has a Missense_Mutation as well as a Frame_Shift_Del and you would like to plot its value, of course you will have a category and color that is different than those of Missense_Mutation or Frame_Shift_Del, that is only natural.
What you can do is to create additional variables like your_gene_Missense_Mutation or your_gene_Frame_Shift_Del each with 0 or 1 to indicate if a particular mutation type is present or not for a particular gene (you can think of it like one-hot-encoding). With a lot of different levels, the resulting heatmap could be difficult to interpret though.
Another simpler alternative would be a using a "universal label" if your gene has multiple mutations. Basically if the cell contains a ;, its value will be changed to something like mixed_mutation_status and instead of multiple "mutation_type_1; mutation_type_2; ... like categories, you will just have one mixed_mutation_status on top of Missense_Mutation, Frame_Shift_Del, ...


Changing the record id in a FASTA file using BioPython

If I use SeqIO.parse(filehandle, 'fasta') to parse a FASTA file, then it will return a SeqRecord object where the id and name are the first word (everything before the first whitespace) of the line beginning with > and the description is the complete line (all not including the initial >). (This behaviour can overruled by providing a custom title2ids function).

So, if I for example have a original.fasta file like:

>Lorem|ipsum>dolor sit amet
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

Then I will get:

>>> from Bio import SeqIO
>>> path = r'C:\path\to\original.fasta'
>>> records = SeqIO.parse(open(path), 'fasta')
>>> record = next(records)
>>> record.id
'Lorem|ipsum>dolor'
>>> record.name
'Lorem|ipsum>dolor'
>>> record.description
'Lorem|ipsum>dolor sit amet'

So, what happens in my case when you do

    if record.id == 'foo':
        record.id = 'bar'

is that the record.id is successfully changed from foo to bar, but that the record.description is not changed and stays foo. That's why, when the FASTA file is printed out, I see the described behaviour of

>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

So, the solution to my problem is to both change the id AND the description:

from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        if record.id == 'foo':
            record.id = 'bar'
            record.description = 'bar'
        SeqIO.write(record, corrected, 'fasta')

with open(corrected_file) as corrected:
    records = SeqIO.parse(corrected, 'fasta')
    for record in records:
        print record.id             # prints bar
        print record.name           # prints bar
        print record.description    # prints bar


How to run 'join' command for multiple files (of 2 types) in a folder

If all of your .tf.txt and .tt.txt files are in the same directory, the following command should work.

for tffile in $(ls *.tf.txt)
do
    prefix=$(basename $tffile .tf.txt)
    join ${prefix}.tf.txt ${prefix}.tt.txt > ${prefix}_final.txt
done


Question about research in bioinformatics for Computer Science student

I think , this kind of research never stop until find way for cure .
I have some research in heart disease,and write this paperECG, you can have some research in improve learning method to find fast and more accurate diabetes , there are lot of database your can search in google and more famous is UCI   .


Where to get '--fldMean' and '--fldSD' for single-end Salmon run

From looking at the salmon documentation here, it is correct that both of these are essential for single end experiments:

  Since the empirical fragment length distribution cannot be estimated from the mappings of single-end reads

From the discussion here it would appear that you can't estimate these parameters at all via the FASTQ files alone.

The fragment length refers to the physical molecule of D/RNA captured in the library prep stage which you have sequenced (in part) on your sequencer.  Thus ideally this sort of data should come from a bioanalyzer trace from whoever performed the library prep or similar QC step, this should ideally show the mean fragment length and standard deviation.


No, downsampling in general has no use as it reduces power. Peak callers such as macs2 will internally scale data to account for sequencing depth differences.
The IgG itself is only useful for peak calling as it somewhat corrects for the unspecificity of the antibody to precipitate any kind of protein rather than the protein of interest. I hear people say, and personally agree with this by experience with my own experiments, that of course IgG pulls down notably fewer DNA that the target antibody, therefore the library requires notably PCR to produce enough material for sequencing. And this then makes it questionable how representative IgG is as a control. I personally prefer chromatin input to account for the intrinsic sequence/GC bias of the underlying genomic region in the sequencing/mapping. To my knowledge there are currently few established tools to really take these controls into account for downstream analysis such as differential testing via interaction terms because the library composition is simply so different from the IPs making normalization difficult.
Your focus imho should be rather to properly normalize your samples to correct for differences is signal/noise ratio, see for example my post over at biostars on how one could do that beyond naive per-million scaling.


How to find hairpin loop structure in a large set of long non coding RNA transcripts

For the first question: "how many of long non coding can make hairpin loop structure?"

The easiest thing to would be to run your sequences through an RNA secondary structure prediction tool. There's a few tools for doing this but the most commonly used are RNAfold from the ViennaRNA package and mfold.

This will give you an output that looks like this, where the parentheses indicate which bases are paired with each other.

GGGCUAUUAGCUCAGUUGGUUAGAGCGCACCCCUGAUAAGGGUGAGGUCGCUGAUUCGAAUUCAGCAUAGCCCA
(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).

To see if there's hairpins, you'll need to annotate the elements within it. Take a look at this question on biostars for an example of how to do this. In short, you can use the forgi library to produce an annotation like this:

(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).
sssssssmmsssshhhhhhhhhssssmssssshhhhhhhsssssmmmmmssssshhhhhhhsssssssssssse

If there's an h in the output annotation, you have a hairpin. The remainder of your question about annotating miRNAs will be left to somebody else.


R packages for data analyses of pooled CRISPR screens

Assuming you mean CRISPR screens targeting many loci (for example using the GeCKO library) there is an R package here.

It uses a linear mixed effect model to compare guide counts before and after a selection step, allowing for multiple guides/gene and multiple replicates.
It can also do the initial read alignment AFAIK and the author should be quite responsive in case you have questions.

Another option could be the python program MAGeCK


Error while loading a big file with Cytoscape

The first thing you can try is to separate the network from the meta information. Assuming you have access to a *nix machine, you can run:

awk -F"," '{print $2,$4}' file.csv > network.csv

That will produce a file like this (based on your example data):

test1 test2
test11 test12
test12 test22

That should radically decrease the file size and will give you the best chance of loading it in Cytoscape. You can then similarly extract the attributes:

awk -F"," '{print $2,$1}' file.csv > source.column1.attributes
awk -F"," '{print $2,$3}' file.csv > source.column3.attributes
awk -F"," '{print $4,$5}' file.csv > target.column5.attributes
awk -F"," '{print $4,$6}' file.csv > target.column6.attributes
awk -F"," '{print $4,$7}' file.csv > target.column7.attributes
awk -F',' '{print $2,$4,$8}' file.csv > edge.column8.attributes
awk -F',' '{print $2,$4,$9}' file.csv > edge.column9.attributes
awk -F',' '{print $2,$4,$10}' file.csv > edge.column10.attributes
awk -F',' '{print $2,$4,$11}' file.csv > edge.column11.attributes
awk -F',' '{print $2,$4,$12}' file.csv > edge.column12.attributes
awk -F',' '{print $2,$4,$13}' file.csv > edge.column13.attributes
awk -F',' '{print $2,$4,$14}' file.csv > edge.column14.attributes

Now, try loading network.csv into Cytoscape and then load each of the node and edge attributes separately.

However, a network with 200,000 nodes and 230,000 edges is just too large to be able to manipulate in any useful way in a GUI. It will just look like a horrible hairball. I really urge you to do some preprocessing and pare the network down to something more manageable. I used to work with a human PPI network of ~13000 nodes and ~15000 edges. While I could load that into Cytoscape, that was already far too big to be able to do anything useful with. I could use the data analysis aspects of Cytoscape, but not really the visualization.


Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?

Bwa-mem may produce chimeric alignment, where different parts of a read are mapped to distinct loci. flagstat counts them as two reads.


Hypergeometic test for the overlap of two set of genomic intervals (e.g. from Chip-seq data)

There is no such thing as a hypergeometric test, at least in statistical textbooks. It's a fisher test based on hypergeometric distribution.

If it is chip-seq for the same target, i.e biological replicates, significance of overlap is not quite meaningful. You get more information by for example checking the correlation of the coverage between your replicates using deeptools or anything else.

If you are comparing between two different targets, for example, chip for two different histone marks, you can use the fisher test or chi-square this way.

Partition the genome in X bp bins, these are all your marbles
Identity those that overlap with peaks from target 1 and those that don't
Do 2. for target2
construct a contingency table using 2 and 3.

So you should end up with something like this (purely made up):

                      no_overlap_target2 overlap_target2
  no_overlap_target1               9799             112
  overlap_target1                    88               1

And if there is overlap, the bottom diagonal with be higher.. less on the off-diagonal, and you can test it with a fisher.test.


As the error tells you: the GROMACS binary has been built without GPU support.
Did you install gromacs with apt-get? Usually these kind of programs are shipped with a basic built, with a lot of features disabled.
If you did built gromacs from source, then cuda was not found during configuration.
I would try building gromacs from source and:

be sure cuda is in your path and check that nvcc and the drivers have a matching version
run cmake and check the output, to see if the video card is recognised
run cmake with -DGMX_GPU=ON

Note that it seems that gromacs doesn't support multiple gpu accelerations( e.g. it is either openCL or cuda, you can't have both)


Drug with as many targets as possible through a system biology graph?

I don't know any python (or R package) to do this), but this is how I would do it based on these matrices:

Do a multidimensional scaling (MDS) plot of the target_similarity matrix.  Explore the visualization:
Find which drug targets the drug closest to the 0,0 point in the 2D space

However this doesn't say anything about the affinity between the drug and the target.
If you wanted to add this then you could weight how much close is the distance between the drug and the target vs how close are they to the center of the MDS (a simple division could work if the distance is scaled between 0 and 1 and averaged).

Note that depending on the original information used to calculate the similarities and the affinities it could be done in other ways


Bioinformatics flowchart for Materials and Methods

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.


awk for motif splicing from an amino acid alignment

I think you are trying to use awk to get characters 3 through 7 (inclusive) of each input sequence. If so, a working version of your approach would be this:

my $splice = `echo "$string" | awk 'BEGIN {FS=\"\"} {x = $n; while (x < $z){ printf \"%s\",\$x; ++x} print ""}' 2>1`;

That, when added to your perl script, produces:

$ foo.pl
VLVIL
VLVIL
VLVIL

That isn't the output you show, but the output you show is just the first 5 letters of each sequence, which doesn't seem to correspond to what your script is doing.

Now, the reason your version didn't work, is that the while was like so:

while (x < $z); print \$ x; ++x}

Assuming the \$ x was a typo, the problem is that the ++x isn't part of the while loop. So that should never even have exited. Seeing as we're in the same time zone, and given the time, I am guessing you're tired and that you miscopied something when trying to give us a reproducible example. In any case, to get the loop right, you need {}:

while (x < $z){print \$x; ++x}

However, if all you want is a substring, a loop is a very inefficient way of doing that. You can instead use awk's substr command:

  substr(s, m[, n ])
             Return  the at most n-character substring of s that begins at position
             m, numbering from 1. If n is omitted, or if n specifies  more  charac‐
             ters than are left in the string, the length of the substring shall be
             limited by the length of the string s.

So, in your script, that would be:

my $splice = `echo "$string" | awk '{print substr(\$0,$n,5)}'

Since you're already using a language that excels in string manipulation, I really don't understand why you would try to use awk instead. It will not be faster or more efficient than the same thing in perl:

print substr($_,$n,5) . "\n" for @array;

By calling awk you are also calling the shell which will be slower by definition, especially if you're echoing thousands of lines. In fact, I tried this with 15000 sequences (just repetitions of your first one) and using either the awk loop I gave above, or the awk+substr, or the native perl substr:

awk loop, average time taken to run the same thing 100 times: 0.0525
awk substr, average time taken to run the same thing 100 times: 0.0522s
perl substr, average time taken to run the same thing 100 times: 0.0464

And you can make it even faster by skipping the array and iterating over the keys directly:

print substr($dna{$_},$n,5) . "\n" for keys(%dna);

That averaged at 0.0376s over 100 runs.


You wrote,
mylist = re.match(H149xcV\_\w+\_\w+\_\w+)

try...
myregex = re.compile(r'H149xcV\_\w+\_\w+\_\w+')
mylist = myregex.match(myfile.id)

or
mylist = re.match(r'H149xcV\_\w+\_\w+\_\w+',myfile.id)

Watch your indentation it should be 4 spaces. Note you are losing the last '_\w+', but don't think it matters.
I think it is better to use,
mylist, _, _ = re.split(r' ', myfile.id)

The problem with your code is list will just return true ... and that is not what you want. If you are using split carefully watch the number of ' ' in your code you can also use r'\s'
You can also set maxsplit as a parameter, which might be useful.

myfile not defined.
Please try,
fasta_file = (sys.argv[0])

Generally your argv command isn't cool and not just [0] versus [1]. I forget the argv module everyone uses .. there's two of them, they are a bit fiddle but are very cool of options. You way in is okay, okayish, sort of okay ... Anyway you appear to be requesting the second input into the commandline, wherea you are probably placing the fasta file in the first position.
Also
print (fasta_file)

Check something happens that you picked up the file.
The other way in is,
mypath = '/User/username/location/fasta.fa'
for myfile in SeqIO.parse(mypath, "fasta"):

where mypath is a OS X style, you're probably using Linux so its... /home/username/pathtofile/fastafile.fa (something like that) or
for myfile in SeqIO.parse('/User/username/location/fasta.fa', "fasta"):

I forgot .... don't ever use list as a variable, by mylist is fine. Again I've not debugged my code.


DESeqDataSetFromTximport all(lengths > 0) is not TRUE

I found this same error as of today using this code:
pad.gene.data <-
  tximport(files = file.list.for.gene.analysis,
           type = "rsem",
           txIn = FALSE,
           txOut = FALSE)

dds.gene.level <-
  DESeqDataSetFromTximport(txi = pad.gene.data,
                           colData = pad.characteristic.data,
                           design = ~ condition)

I had done this same differential gene expression analysis using transcript counts aggregated to the gene level, on a mac about 5 weeks ago, but as of this past week switched to a windows and re-installed R 3.6.3 and whatever the current version of tximport and DESeq2 is.
I was able to replicate the results I had gotten on the mac 5 weeks ago removing all rows where all gene lengths are greater than 0:
pad.gene.data$abundance <-
  pad.gene.data$abundance[apply(pad.gene.data$length,
                                1,
                                function(row) all(row !=0 )),]

pad.gene.data$counts <-
  pad.gene.data$counts[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

pad.gene.data$length <-
  pad.gene.data$length[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

So first removing any gene where there is a length count of 0 from the abundance and counts matrix, then removing the genes from the length matrix so that they are all the same size. It was able to run that way.
I am definitely skeptical of why I would need to do that as it doesn't seem that it was advertised, if this is a new version, that we need to preprocess the data beforehand. I hope someone who might know more of the code can comment on this.


how to estimate phylogeny without DNA sequences?

Sequence data is basically always better and at this point in history it is de rigueur for any serious study, but originally phylogenies were all based on (ideally discrete) phenotypic characters, so there's no reason why it can't work in principle. As @Pallie mentions in a comment, the clustering methods in e.g. NJ or UPGMA will work on any set of distance measures, and they were in fact originally developed on phenotypic characters rather than sequence data.
I'd recommend taking a look at chapter 24 of Joe Felsenstein's book on phylogenies, which is digitized here (yes the link looks sketchy, sorry). There is a section on p. 425 "Inferring phylogenies and correlations", which discusses how you can in principle use quantitative characters to infer phylogenies. The theory and math for why this is supposed to work is earlier in the chapter.
In many cases people prefer to code the characters as discrete because it makes it easier to work with. But there are also methods that work with the raw quantitative characters.
In R, you should be able to do something like this:
# some simulated data, species are rows and characters are columns
> data
      char1 char2 char3
spec1     0     3     6
spec2     1     4     7
spec3     2     5     8

> library(ape)
> tr = nj(dist(data))
> plot(tr)


data storage requirements for RNA-seq and WGS

The breakdown would be:

1.5-2 GB per sample for fastq files.
~3GB per BAM file, depending on the aligner

So ~5GB per BAM file, which is about 50GB total. That can vary by maybe 20%, depending on how compressible the results are (I rearrange our fastq files to make them more compressible, so I'm probably underestimating what others would see).

Obviously using Salmon/Kallisto would halve the storage requirements. Using CRAM would also decrease things, though obviously fewer tools support it.


Unable to install bedtools on windows 10 ubuntu

My personal advice would be to avoid compiling tools like this if you can avoid it, especially where you are using a "barebones" setup like the one you get with WSL Ubunutu as you'll almost always run into dependency issues.

The linux version of conda installs just fine in WSL and bioconda hae pretty much any bioinformatics package you might need, already compiled.


Remove variants that do not map to human genome

Answer I received on Biostars that appears to work as expected:

bcftools norm -Ob --check-ref wx -f <reference_fasta> <vcf.gz> -o <output_name>


Looking at this resource, these data are microarray data rather than RNAseq data. Thus the "signal" in question is the intensity of florescent light being emitted from the the probe associated with your gene.  The relationship between intensity and quantity of RNA is different for every probe because each probe will have a different affinity for complementary fragments dependent on the sequence in question (not least the melting temperature).

Affy intensities are usually presented as log values going up to around 12, but the raw number is next to meaningless and usually only means something as a normalized (e.g. RMA) log fold change.

I've never come across RefExA before, and it appears to be undocumented and unmaintained. I wouldn't rely on it for usable data if I were you. A more modern database of expression data is the EBI Expression Atlas. It presents data from a wide range of different studies, processed and normalised on a unified pipeline and expressed as TPMs.


Total reads aligning to each reference within a bam file

The quick way to get the number of alignments on each reference is

samtools idxstats my_bam.bam

Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us this as generally I'm after a rough approximation, rather than an accurate number.

In theory, only one alignment for each read should be marked as primary, so the following should give you what you need quickly and at low memory usage:

samtools view -bF 2304 my_bam.bam > primary_only.bam
samtools index primary_only.bam
samtools idxstats primary_only.bam


How to extract values from second file on the basis common first column?

I like to use the join command for this. For your example above, you can simply use the following:

$ join file2 file1
ENST00000000412 0.779142 M6PR
ENST00000001008 0.738143 FKBP4
ENST00000002501 0.715315 DBNDD1
ENST00000002596 0.713664 HS3ST1

This assumes your files are sorted. If they aren't sorted, you can use input redirection if your shell supports it (bash and other modern shells do):

join <(sort -k 1b,1 file2) <(sort -k 1b,1 file1)

If your shell doesn't support <(), sort the files first:

sort file1 > file1.sorted; sort file2 > file2.sorted
join <(sort -k 1b,1 file2.sorted) <(sort -k 1b,1 file1.sorted)


Trying to divide number in fourth TSV column by a certain number according to its corresponding identifier in the first column

Python dictionaries are indeed a good tool for that.

Here is a script that takes 2 arguments (the file containing the counts and the file containing the genome length info), and prints the input data, with an added column containing coverage:

#!/usr/bin/env python3
import sys

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

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

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

sys.exit(0)

Testing it:

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


How to retrieve logical expressions (KO based) for reactions from KEGG?

Some information seems to be available on the module level utilising the API.

Reaction R00352 points to the module M00173. While some information is in the Definition, it's not trivial to link it to the particular reaction. This linkage is however easier accessible via API:


Counting the characters in a row using python


Why do NEB adapters have non-complementary sequence?

The original point of the non-complementary sequence is to allow formation of a Y-adapter after USER treatment, since that's what's needed for the Illumina chemistry to work. In reality the bottom strand may not be there in step 3, since USER has already digested one strand (creating the hairpin in step 2) and size selection has been performed. Regardless, it'll get removed once cluster formation is performed.


How do I prevent the FeatureHeatmap function from the Seurat package, from sorting my data groups in alphabetical order when plotting data?

I had this same issue. I reordered the levels and that worked for most functions, but would switch back with the FeatureHeatmap. Turns out there is a line in the function that changes it to a character.

I ended up making a "new" function and just changed that one line from:

data.plot$ident <- as.character(x = object@ident)

to:

data.plot$ident <- object@ident

Once that was fixed, the order that I assigned worked out just fine.


How can I get a list of papers where a particular mutation has been mentioned?

You could use the pubmed API to query it directly.  Typically I use Biopython for my pubmed queries.  There are some good examples in the cookbook and on other SO posts.

This example from the cookbook link above has the guts of what you need.  You will just have to iterate over your list.  Depending on the number of queries you want to do you might have issues getting throttled by NCBI

In this example, we will query PubMed for all articles having to do with orchids (see section 2.3 for our motivation). We first check how many of such articles there are:

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="orchid")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]: ...
        if row["DbName"]=="pubmed": ...
            print(row["Count"]) 463

Now we use the Bio.Entrez.efetch function to download the PubMed IDs of these 463 articles:

>>> from Bio import Entrez
>>> handle = Entrez.esearch(db="pubmed", term="orchid", retmax=463)
>>> record = Entrez.read(handle)
>>> handle.close()
>>> idlist = record["IdList"]

This returns a Python list containing all of the PubMed IDs of articles related to orchids


In Deseq2, it depends on contrast eg. if your contrast <- c("Condition", firstC, SecondC) is then -ve is downregulated in FirstC and same for +ve. It is always better to take look at count data for the confirmation (responders versus non-responders samples).
[for more detail you can refer here][1]
[1]https://github.com/amarinderthind/RNA-seq-tutorial-for-gene-differential-expression-analysis/blob/master/RNA-Seq-DGE.R


This sort of procedure is done in tool papers to benchmark results compared to different aligners/mappers. See the salmon or kallisto papers for examples. Using the intersect of DE results from different mappers will end up yielding somewhat conservative results, not unlike using different differential expression packages and taking the intersect of their outputs. In practice you could just take a more conservative p-value cut-off and save the hassle of using multiple tools, the results would end up being similar.


Without giving away too much for your homework check out the sample() function in R. Think about how you can pass DNA bases to this function to generate random sequences.


Every version of SBOL allows you to specify prefixes/suffixes as annotations on a sequence.

In SBOL 3, the construct would be a Component of type DNA (https://identifiers.org/SBO:0000251), and the linker a SequenceFeature on that Component. The role of the sequence feature might depend on the particular assembly you are describing, and whether you're talking before or after assembly (in those cases where the process changes the sequence), but a good one in many cases would be Nucleotide Binding Site (https://identifiers.org/SO:0001655)
In SBOL 2, the representation is the same, except the construct would be a ComponentDefinition and the linker a SequenceAnnotation.
In SBOL 1, the representation is still the same, except the construct would be a DnaComponent and the linker an SBOL 1 SequenceAnnotation (which has different fields than in SBOL 2).

If the assembly method does not change the sequence, then all three can represent this pretty much equally well.  If you want to represent a sequence that is modified by assembly or the plasmid carriers from which parts are extracted, then SBOL 3 is the best because it cleanly represents deletions during assembly via the sourceLocation field.


While I haven't found a way to limit the results to the canonical transcript only, you can get a list of genes, transcripts and their CDS lengths using Ensemble's BioMart. I have already set it up for you, you can see the results, and modify them, here (click on the "Results" link if you don't see them).

Essentially, you just need to go to BioMart, and

select "Ensembl Genes 96" (the number will change if the version changes) as the database and "uman Genes" as the dataset.
Click on "Filters", and set Gene type to coding and Transcript type to protein_coding.
From "Attributes", select whatever you want to see. The "CDS Length" is under "Structures".


Why is my STAR reference genome indexing aborting on my GNU/Linux server but not on my Mac OS X laptop?

A std::bad_alloc error tends to mean that you've run out of memory. It's not unusual for head nodes to be fairly limited in capacity, so presumably that's the issue.


removing mitochondrial read and unassembled "random" from multiple bam files

Your awk command is inside double quotes ("") which means that all of the awk field variables (e.g. $3) are being interpreted by the shell. Since $3 is not set in the shell, your awk is trying to compare an empty string and that's why you get the error:

awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error

Note how the error message shows if( != chrM and not if($3 != chrM. Since you need to use double quotes here because you need both an outer (for parallel) and an inner (for awk) quote, and single quotes can't be nested, you will have to escape the awk variables:

ls *.bam | parallel "samtools view -h {} |
    awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb {.}.bam"


Check if a set of genes fall more on a specific cluster than other clusters

A hypergeometric test is related to the fisher test, in a table with one  dimension being some factors and the other one some other factors it is tested if they are independent. In R you con do this:

data <- matrix(c(0, 1, 2, 5), ncol = 2, nro = 2)
fisher.test(data)

Or you can use the specific function phyper, to see the relation between them you can read this question and answer


Standard for when to collapse phylogenetically uninformative nodes

If 1E6 is genetic distance, this means one nucleotide difference per million base pairs, which is of course incredibly small and bootstrapping at this level doesn't really work. So what I suspect is happening is that the bootstrap is so low it is not present on the tree. You can assign a cut-off bootstrap but anything below 75% is not interpretable in any case. What might be happening is you are producing the bootstrap "cladogram" which is showing a default branch length, but when <1E6 genetic distance it is showing a trifurcation. This simply means the bootstrap of the consensus tree is so low it is collapsing the branch and this is a function of the very small genetic distance at this node.

Looking at your treefile I don't see any collapsed clades.
(Group_9:0.00000100000050002909,(Group_50:0.00507659947127924935, (.

The clade above has not been colapsed it is just a very small distance in a paraphyletic group. It looks fine to me, i.e. Group_9 is 1E6. It could look 'collapsed' simply because it is so small it is difficult to see. I hope thats ok?


How to compare transcriptomic profiles of two cell types (single cell RNA-seq)?

You are looking for a way to compare expression profiles between cell types, and your data is counts of genes expressed per cell. Your issue is that the data is highly dimensional. It has as many dimensions as you have genes. To be able find useful information, you must reduce these dimensions using a statistical technique such as Principal COmponents Analysis (PCA), t-Distributed Stochastic Neighbor Embedding (t-SNE), or Multidimensional Scaling (MDS). The result of these analyses is a reduction of your thousands of dimensions to vectors (eigenvectors in PCA) of variance which you can plot against each other in a scatter plot.

I have been working with Single Cell RNA-seq data and have been exploring the difference in gene expression in a certain type of cells between timepoints. Seurat (R package) is a very powerful tool for scRNA-Seq data analysis.

Here is a a tSNE of my data. I have added the timepoint labels, but you can see that after performing the dimensionality reduction the two different timepoints cluster by their gene expression profiles. I can then use Seurat to see which specific genes are differentially expressed between the clusters and a whole bunch of other fun things.


Why do some assemblers require an odd-length kmer for the construction of de Bruijn graphs?

From the manual of Velvet:

  it must be an odd number, to avoid palindromes. If you put in an even
  number, Velvet will just decrement it and proceed.

the palindromes in biology are defined as reverse complementary sequences. The problem of palindromes is explained in this review:

  Palindromes induce paths that fold back on themselves. At least one
  assembler avoids these elegantly; Velvet requires K, the length of a
  K-mer, to be odd. An odd-size K-mer cannot match its reverse
  complement.

It is possible to construct graph with palindromes, but then the interpretation will be harder. Allowing only graphs of odd k-mers is just an elegant way to avoid writing a code for interpretation of a more complicated graph.


How do I infer someone's haplotype of a gene from data on several of the gene's SNPs?

It depends on the scale of the data (i.e. number of individuals, number of SNPs). Generally this can be done by matching up individual SNP genotypes to haplotype tag SNPs. The specific combination of genotyped SNPs is used to determine the underlying haplotype, based on tag SNP frequencies for each haplotype.
More details of the process of identifying tag SNPs can be found here:
https://www.broadinstitute.org/haploview/tagger


How is the GT field in a VCF file defined?

0 and 1 is just a way of coding the reference (0) and alternate (1) allele. These could be A/G, C/T etc. It's just a simplified way of expressing the different alleles.
0/1 and 1/0 functionally mean the same thing (that the individual is a heterozygote) - since the genotype is unphased, the alleles aren't ordered. The / symbol tells you the genotype is unphased. However, in practice, a heterozygous genotype is always written as 0/1 as a matter of convention. 0/0 is referred to as homozygous reference, and 1/1 as homozygous alternative.
0|1 and 1|0 do mean something different, since the pipe symbol | tell us the order of the alleles matters.
Occasionally you will have multi allelic positions, where you have more than one alternate allele, and therefore the field will look like:
Ref   Alt
A   [C,GT]

The extra alt allele is given the number 2, so in this case if an individual had the genotype CT then their genotype code would be 0/2.
Full detail is given in the official VCF specification page.


How to calculate Gene Ontology terms in python

You can calculate gene ontology similarity with the GOSemSim R package (paper)(disclaimer I'm a contributor to this package). You have several similarity scores implemented, some of them are across the three subontologies.

In python you have the GOATOOLS (paper). I think there was another Python package but now I don't remember its name.


Extracting expression data from GSE dataset downloaded from GEO

According to the manual, all you need to do is:

library('GEOquery')
gseGSE16146 <- getGEO('GSE16146', GSEMatrix=FALSE)

As explanation, getGEO() outputs by default to GSEMatrix=TRUE and returns a list of ExpressionSet objects. You should get what you were looking for with:

Table(GSMList(gseGSE16146)[[1]])[1:5,]

The manual has also a paragraph about this:

  The GSE has  a  metadata section, just like the other classes.
  However, it doesn’t have a GEODataTable.  Instead, it contains two
  lists, accessible using GPLList and GSMList, that are each lists of
  GPL and GSM objects.

I tested the code I posted on R 3.3.2 and I get access to the Data Table just fine.

update:

Answering how to extract the actual expression data (not what asked in the original question):

library(GEOquery)
data = getGEO("GSE16146")
datExpr = exprs (data[[1]])


hgsql not found when running a script despite being accessible on command-line

It is possible that the shell being used in the process is different from your standard command line shell. I am not familiar with the tool, and it's hard to say more without knowing your OSX version, but it looks like it is using zsh. I think that zsh only fairly recently became the standard OSX shell. I don't know the tool or the shell very well, but it looks like you might be able to explicitly set up hgsql in your PATH in your .zshrc.

It's also possible that the liftover script was not built to run via zsh. I notice that it is trying to run csh scripts as part of the workflow. Then again the localhost may insulate it.

So I think that my first step would be to sort out which shell I am using, whether it works with the liftover pipeline, and how it's running on the localhost. I notice that the ssh call running the csh script on the localhost is the step that actually fails. What shell is executing the csh script? Does that shell know about hgsql? It could be that the way that the error bubbles up it is hiding which shell it is actually using.

I'd recommend trying that ssh call in isolation from the rest of the pipeline, and diagnosing the error that way.

Hope that helps.


I think what makes this especially hard is the way your data is formatted.
The problem would be way easier if you had a matrix of presence / absence instead. Each row would be a gene, each column a sample, and the values would be 0 (absent) or 1 (present).
Awk becomes a bit inconvenient for this task, so I'll use python instead. To get such a matrix from your file, you could do the following:
import numpy as np
import pandas as pd

# Get total number of genes / samples
all_samples = set()
all_genes = set()
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    all_genes.add(line[0])
    all_samples.update(line[2:])
all_samples = list(all_samples)
all_genes = list(all_genes)

# Initialize presence matrix with correct dimensions
pres_mat = np.zeros((len(all_genes), len(all_samples)), dtype=int)

# Map each sample/gene to a row/col
row_map = {gene: i for i, gene in enumerate(all_genes)}
col_map = {sample: i for i, sample in enumerate(all_samples)}

# Fill matrix incrementally
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    curr_gene = line[0]
    row = row_map[curr_gene]
    for sample in line[2:]:
      col = col_map[sample]
      pres_mat[row, col] = 1

# Add names and write to file
df = pd.DataFrame(pres_mat, columns=all_samples)
df.insert(0, 'gene', all_genes)
df.to_csv('gene_presence.tsv, sep='\t')

You should get something that looks like:
       gene  IZ_BH_5.5  A8.S8  IZ_BH_2.1  A9.S9*  IZ_BH_5.3  A9.S9  ...  IZ_BH_3.11  IZ_BH_4.7  IZ_BH_4.9  IZ_BH_2.8  IZ_BH_3.5  IZ_BH_2.6  IZ_BH_4.3
0       Lyn          1      1          1       1          1      1  ...           1          1          0          0          1          1          1
1      Nfib          1      1          1       1          1      1  ...           1          1          1          1          1          1          0
2  Cdk5rap2          1      1          1       1          1      1  ...           1          0          1          1          1          1          1

With this format, it is now easy to compute the most frequent samples. E.g. df.sum(axis=0) to get the number of occurence of each sample for these 3 genes. It would also allow you to do different things, such as computing the pairwise similarity between samples.


Samtools/bcftools calling indels in noisy reads

It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021.
What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the pileup calls for the surrounding bases may help to diagnose this.
My general advice to people doing variant calling with nanopore is to be suspicious of deletions, especially if they include homopolymers. That usually means trying to confirm any important variants with another technology.
Newer, more accurate, base calling models are being created all the time by ONT, so if you still have the raw FAST5 data, it might be worth it to run through the most recent base caller.


How to find the position-specific weight matrix (PSWM) for a list of genes having different sizes?

If your sequences are different lengths, this implies a gapped alignment with insertions/deletions. This is not possible with a regular PSSM, so you have to create a separate PSSM for each group of homologous sequences of the same sequence length.

The alternative is to use a different approach that allows gaps like profile HMMs (as implemented in HMMER3 or HH-Suite3).


Comparison of gene set enrichment statistics

I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply.

I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:

library("limma")
set.seed(123)

# Create some genes and samples
nGenes <- 40
nSamples <- 10
expr <- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes <- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) <- genes[seq_len(nGenes)]
colnames(expr) <- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values
expr2 <- expr+0.5
rownames(expr2) <- genes[seq(41, to = 80)]
expr <- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] <- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")

# Create a design matrix and
design <- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)),
                 ncol = 2, nrow = 10)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
tt <- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet <- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
##
## $correlation
## [1] -0.02036983
##
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
##
## $correlation
## [1] 0.3024643
##
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
##
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392

In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance.

As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:

               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97

As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.

  How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value?

I am unsure how to interpret that comparison. It would be better to explore this  with real data (for example from ALL library), and contact the authors, either via the support site in Bioconductor or through personal communications.

  Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?

As you might have noticed I used inter.gene.cor = NA in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower.

Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (cameraPR) because it doesn't have the same results as in camera and it doesn't allow to estimate the intercorrelation.


Counting a specific consecutive character with its occurrence position and length

I couldn't come up with a one liner but have a two liner instead...

To get the start positions of the N repeat:

grep N mySeq.fasta | grep -b -o N | awk -F: '{print $1}' |awk '$1!=p+1{print $1}{p=$1}'
12
38
56

To get the length of the repeats:

grep N mySeq.fasta | awk -F '[^N]+' '{for (i=1; i<=NF; i++) if ($i != "") print length($i)}'
14
6
9

With thanks to stack overflow.


How to download FASTA sequences from NCBI using the terminal?

Alternatively, you can use the NCBI Entrez Direct UNIX E-utilities

Basically, you have to download the install file here: https://www.ncbi.nlm.nih.gov/books/NBK179288/bin/install-edirect.sh

In the terminal, install it using: source ./install-edirect.sh

Then, you can download your sequence by doing:

esearch -db nucleotide -query "NC_030850.1" | efetch -format fasta > NC_030850.1.fasta

And you should find your fasta sequence downloaded.

As you have several sequences to download, I think it will be quite easy to add this command into a little bash script to process all of them.

More readings about that: http://bioinformatics.cvr.ac.uk/blog/ncbi-entrez-direct-unix-e-utilities/ and https://www.ncbi.nlm.nih.gov/books/NBK179288/


Plink 1.9 --merge-list changing order of fam file

Thanks Christopher Chang for the great answer via the plink2 google group! (See here)

His answer was as follows:

"Yes, plink 1.9 did change the default merged sample order.  However, you can request plink 1.07’s behavior by adding “--indiv-sort 0” during the merge.  (You can also use --indiv-sort + --make-bed at any time to change the sample order to whatever you want.)

As for why the default was changed, take a look at the new default order: it will probably make a lot of sense.  It’s “natural-sorted” which tends to make it easy for humans to manually search through regardless of what mix of letters and numbers are used in the IDs."

See also the Plink 1.9 documentation for --indiv-sort which makes the point that this behavior is default only for merges--in other commands Plink 1.9 uses the same ordering as in Plink 1.07, which is the order the samples were in.


can I download DESeq2 in R 3.6.3 in Linux MInt?

Bioconductor is a package management tool for R that allows you to install the appropriate package versions for your R version.
For R 3.6.3 you'll need Bioconductor 3.10:
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.10")

Then you install  DESeq2:
BiocManager::install("DESeq2")


How to remove sequence reads contains more than 2 X from multifasta file?

A Biopytnon solution that works for .gz files too :-). The script make a script filter_faa.py:

#!/usr/bin/env python3

import sys
from Bio import SeqIO
from mimetypes import guess_type
import gzip

filename = sys.argv[1]

encoding = guess_type(filename)[1]
_open = partial(gzip.open, mode='rt') if encoding == 'gzip' else open

with _open(filename) as ffile:
    for seq_record in SeqIO.parse(ffile, 'fasta'):
        if seq_record.seq.count('X') < 3:
            print(seq_record.format("fasta"))

Then usage is

python3 filter_faa.py my.faa > my_filtered.faa

The my_filtered.faa then will contain

>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF

>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

The filtering condition is hard-coded. You could potentially parametrise it if you feel like later you might want to filter by a different number of unknown aa.


Getting all NNI trees of a parsimony tree

As I've studies the 'TreeConstruction' module of Biopython, I came to know that '_get_neighbors' method is the one that is used to get all neighbor trees of the given tree. This is all I've concluded from my findings. If anyone wants to add more, that'll be positively appreciated.


I found a post useful for this topic.

It explains the difference of coverage and depth.  It also has a useful explanation on how to calculate coverage and depth.

Here is a copy of what the link says just encase the post is removed:

  Depth of coverage

  How strong is a genome "covered" by sequenced fragments (short reads)?

  Per-base coverage is the average number of times a base of a genome is
  sequenced. The coverage depth of a genome is calculated as the number
  of bases of all short reads that match a genome divided by the length
  of this genome. It is often expressed as 1X, 2X, 3X,... (1, 2, or, 3
  times coverage).

  Breadth of coverage

  How much of a genome is "covered" by short reads? Are there regions that are not covered, even not by a single read?

  Breadth of coverage is the percentage of bases of a reference genome
  that are covered with a certain depth. For example: 90% of a genome is
  covered at 1X depth; and still 70% is covered at 5X depth.

SAMtools: get breadth of coverage

Get breadth of coverage of a reference genome.

How many bases of my reference genome are covered by my short reads sample at a level of 5X or higher?

 # set shell variables to your file-names

 REF_GENOME_FILE=my_ref_genome.fna
 WGS_SAMPLE=my_short_read_sample.fastq
 MIN_COVERAGE_DEPTH=5

Map short reads against reference genome (bowtie2)

# create bowtie2 index database (database name: refgenome)

bowtie2-build ${REF_GENOME_FILE} refgenome
ls
refgenome.1.bt2
refgenome.2.bt2
refgenome.3.bt2
refgenome.4.bt2
refgenome.rev.1.bt2
refgenome.rev.2.bt2

# map reads and sort bam file
bowtie2 -x refgenome --no-unal -U ${WGS_SAMPLE} -S - -p 12 | \
samtools view -bS - | \
samtools sort -m 5G - mapping_result_sorted.bam

-U ${WGS_SAMPLE} consider all reads as unpaired reads
-S -             bowtie output in SAM format, write to standard out
-x refgenome     use reference genome bowtie2-database
--no-unal        suppress SAM records for reads that failed to align
-p 12            use 12 processors

Get breadth of coverage (SAMtools)

# create samtools index
samtools index mapping_result_sorted.bam

# get total number of bases covered at MIN_COVERAGE_DEPTH or higher
samtools mpileup mapping_result_sorted.bam | awk -v
X="${MIN_COVERAGE_DEPTH}" '$4>=X' | wc -l
32876

# get length of reference genome
bowtie2-inspect -s refgenome | awk '{ FS = "\t" } ; BEGIN{L=0}; {L=L+$3};
END{print L}'
45678

Result: breadth of reference genome coverage
total number of covered bases: 32876  (with >= 5X coverage depth)
Depth of coverage (average per-base coverage): 0.719 X (32876 ÷ 45678) (total number of covered bases divided by reference genome length)
percent: 71.9%  (0.719 × 100) (coverage breadth in percent)


How are snakemake's --cluster and --drmaa options implemented?

I'd always kind of wondered how this worked too, so I took this as an excuse to look into the snakemake code. At the end of the day this becomes a question of (1) how are jobs actually submitted and (2) how is it determined if jobs are done (and then whether they failed)?

For DRMAA, python has a module (appropriately named "drmaa") that wraps around the libdrmaa library that comes with most schedulers. This is a very popular route, for example the Galaxy project uses this for dealing with most clusters (e.g., I use it to connect our internal Galaxy instance to our slurm cluster). The huge benefit here is that DRMAA does magic to allow you to simply submit commands to your cluster without needing to know whether you should execute qsub or srun or something else. Further, it then provides methods to simply query whether a job is running or not and what its exit status was.

Using the --cluster command requires a LOT more magic on the side of snakemake. At the end of the day it creates a shell script that it then submits using the command you provided. Importantly, it includes some secret files in that script that it can then watch for (ever noticed that .snakemake directory where you execute a script? That appears to be among the things it's used for.). These are named {jobid}.jobfinished and {jobid}.jobfailed and one of them will get touched, depending on the exit status of your command/script. Once one of those is there, then snakemake can move on in its DAG (or not, if there's a failure). This is obviously a LOT more to keep track of and it's then not possible for snakemake to cancel running jobs, which is something one can do easily with DRMAA.


Translating a genome sequence into possible frames

The 6-frame translation  results are below you will need to do the rest. It is translated 5' to 3' but appears to contain two splice sites. You will need to use blast to resolve the start/stop of the splicing.

All 6-codon translation is available the Expasy server here. You will need to blast to resolve the splicing.

5'3' Frame 1
KAIIQFKTWVSRGVRTWCWVGAGRGKGRRGCVRFSWNR-LGKTRANLRTQP-LPCRGRPRVSSRPRCGGTARSGALGGSRTAATLGGVAPGEVSQRGLWRSCISVLSAKGGAPIII-HPPCIYGGVLKPAAELATPAGERKKKSWQKECWTGAVLGVGTGAVEREGWEGLRCRRG-ESG-GASGNSRSGGAPARSGVHAAASAQAPLAFSFRWRKTLCLGFWQIVFLTATSRGFLRAPGPISIPLPLRGRLPGFALRAPGGAGARRAPSRWFTKCVSEIAGDCPKGVKGCSLYSPTKTTQPL-GIALQGERFGTVARTARCARFPHQDPFNSRLPPALCAPLQQPPATMPLNVSFTNRNYDLDYDSVQPYFYCDEEENFYQQQQQSELQPPAPSEDIWKKFELLPTPPLSPSRRSGLCSPSYVAVTPFSLRGD
5'3' Frame 2
KQ-YNLKPGSLEVLGRGVG-AQAGEKGGEDVSDSPGIVDLEKPGRISAPSPDSPAAAALGCPRARDAEELRGAGLWAVPEQLLPLVGWLRGRYRSGVSGAVASPY-VRREVPLLLFDTPLVFMEGC-SPRLSSPLQPARERRKAGKRSVGRGRYWGWGRGRWRGKVGRGCGAGGGRRAARARVGTAAAEGPRRGAGFTQPLAPRRLSPSPSGGAKLCALDFGKLFSSPPPPAAS-GRQGRFRFLCRCGADSRALRSGLPGERGLGGHQAAGSLSASPR-QGTVQRG-KGAPFIPPPRPPSRFRG-LCKGRGSGLWRALRAAPGFRTKTPLTQDCLPLCVPRSSSLPRRCPSTLASPTGTMTSTTTRCSRISTATRRRTSTSSSSRASCSPRRPARISGRNSSCCPPRPCPLAAAPGSARPPTLRSHPSPFGET
5'3' Frame 3
SNNTI-NLGL-RC-DVVLGRRRQGKREARMCPILLESLTWKNQGESPHPALTPLPRPPSGVLAPEMRRNCEERGSGRFQNSCYPWWGGSGGGIAAGSLAQLHLRIECEGRCPYYYLTPPLYLWRGVKARG-ARHSSRREKEEKLAKGVLDGGGTGGGDGGGGEGRLGGAAVPAGVGERLGREWEQPQRRGPGAERGSRSR-RPGASRLLLQVAQNFVPWILANCFPHRHLPRLLKGARADFDSSAAAGPTPGLCAPGSRGSGGSAGTKPLVH-VRLRDSRGLSKGGERVLPLFPHQDHPAALGDSSARGEVRDCGAHCALRQVSAPRPL-LKTASRFVCPAPAASRDDAPQR-LHQQEL-PRLRLGAAVFLLRRGGELLPAAAAERAAAPGAQRGYLEEIRAAAHPAPVP-PPLRALLALLRCGHTLLPSGR
3'5' Frame 1
CLPEGRRV-PQRRRASRARSGG-GTGAGWAAARISSRYPRWAPGAAARSAAAAGRSSPPRRSRNTAAPSRSRGHSSCW-S-R-GASSREAAGAGHTKREAVLS-RGLGAETWRSAQCAPQSRTSPLAELSPKAAGWSWWGNKGSTLSPPLDSPLLSRRRT--TSGLVPAEPPLPREPGAQSPGVGPAAAEESKSALAPLRSRGRWR-GKQFAKIQGTKFCAT-RRRREAPGR-RLREPRSAPGPLRCGCSHSRPSRSPTPAGTAAPPNLPSPPPPSPPPVPPPSNTPFASFSSFSRRLEWRAQPRALTPLHKYKGGVK---GHLPSHSIRRCNCARDPAAIPPPEPPHQG-QLFWNRPEPRSSQFLRISGARTPEGGRGRGVRAGCGDSPWFFQVNDSRRIGHILASLFPCLRLPNTTS-HL-RPRF-IVLLL
3'5' Frame 2
VSPKGEGCDRNVGGRAEPGAAARGQGRGGQQLEFLPDILAGRRGLQLALLLLLVEVLLLVAVEIRLHRVVVEVIVPVGEANVEGHRRGRLLERGTQSGRQS-VKGVLVRKPGAARSARHSPEPLPLQSYPLKRLGGLGGGIKGAPFHPLWTVPCYLGDALSEPAAWCPPSPRSPGSPERKARESAPQRQRNRNRPWRP-EAAGGGGEENNLPKSKAQSFAPPEGEGERRLGASGCVNPAPRRGPSAAAVPTRALAALLPPPAPQPLPTFPLHRPRPHPQYRPRPTLLLPAFLLSLAGWSGELSRGL-HPSINTRGVSNNNRGTSLRTQYGDATAPETPLRYLPRSHPTKGSSCSGTAQSPAPRSSSASRARGHPRAAAAGESGLGAEIRPGFSKSTIPGESDTSSPPFSPACAYPTPRPNTSRDPGFKLYYCF
3'5' Frame 3
SPRREKGVTAT-EGEQSPERRLGDRGGVGSSSNFFQISSLGAGGCSSLCCCCW-KFSSSSQ-KYGCTES-SRS-FLLVKLTLRGIVAGGCWSGAHKAGGSLELKGSWCGNLAQRAVRATVPNLSPCRAIP-SGWVVLVGE-REHPFTPFGQSPAISETHLVNQRLGARRAPAPPGARSAKPGSRPRSGRGIEIGPGALKKPREVAVRKTICQNPRHKVLRHLKEKARGAWALAAA-TPLRAGAPPLRLFPLAP-PLSYPRRHRSPSQPSLSTAPVPTPSTAPVQHSFCQLFFFLSPAGVASSAAGFNTPP-IQGGCQIIIGAPPFALNTEMQLRQRPRCDTSPGATPPRVAAVLEPPRAPLLAVPPHLGREDTRGRPRQGSQGWVRRFALVFPSQRFQENRTHPRLPFPLPAPTQHHVLTPLETQVLNCIIA


The paf format does not have the CIGAR string and base quality information that are typically present in a BAM/SAM file. While it is in principle possible to write a paf2sam converter, it's going to be very difficult to regenerate the CIGAR string with all the indel information (and maybe it's probably impossible to do so) that you'll need in a BAM/SAM file. If your motivation is to get coverage information, it might make more sense to convert the paf file into a bed file, before using a tool like bedtools to convert it into a coverage track.


Does read sequence identifier need to be unique for proper alignment?

What tool did you use to replace the read names? Can you modify that to split up the reads into separate files or encode the cell barcode in the read description instead of the read name?

Typically is it assumed by many tools that the FASTQ read names are unique. you could have issues with BWA and split alignments but I don't know if that pertains to your data.  I think the best you can do is the sanity check the results at each step of the process.

In the fastq format (https://en.wikipedia.org/wiki/FASTQ_format) on the first line denoted by the '@' there can be two values separated by a space.  The first is the name and the second is the description.

So an example would be:

@Original_read_name CGAGGCTG_TATCCTCT_AGTCCA
ATCGATCATGATGACA
+
@&%*@(%@%*($!##!


As @Michael said PyMOL can do this. PyMOL is not only an app, but also a python package —installed via conda not via the regular download (you can have both).
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(infile,'myprotein')
     pymol.cmd.save(infile.replace('.cif', '.pdb'), selection='myprotein')

This runs the parallelisable instance, so you can have that context manager block running on how many threads or processes you fancy.


Specific cell type identification in Single Cell Sequencing

I do not know of such software.

However, I believe this effort is a bit misdirected. The purpose of single-cell sequencing is to get a better understanding of cells; their heterogeneity and functional diversity or developmental / biological processes such as differentiation, using a higher "resolution" method. In other words, if we had the methods you are asking for, there would be no need for the single-cell experiment.

Running your data through some pipeline utilizing previous knowledge has the danger of forcing ideas on the data, rather than seeing what the data tells you: it would be better to try and understand and explain the biologically relevant heterogeneity and diversity in your cells, together with a critical comparison of cell type characteristics (expression of genes, pathways) with knowledge described in literature.
This question is also very similar to the one you previously asked and was answered.

Also good to recognize that your approach has a lot of implicit assumptions, such as: cells can be clearly categorized as +/- expression for each gene, or it ignores the systems level (networks and pathways), thus again forcing concepts on the data.

This said, one method would be to construct a table listing cell type markers as described in the linked answer above, then write a script that determines a cutoff value for expression of genes in your data (see this), then ranks cell types for each cell. For example, you can measure the number or proportion of uniquely genes expressed (although that's very simplistic).

Another option, if you are familiar with machine learning, is to train a classifier on an annotated dataset, and then use that on new data.

Also see a convenience function below, which requires a table of marker genes for each cell cluster (i.e. (the Seurat::FindAllMarkers() output), and a reference df listing genes and corresponding cell types in its HGNC_symbol and Cell_type columns, and returns the table with listing clusters, their marker genes and corresponding cell types.

getCelltypes <- function(markers, reference) {
  marker.celltype <- markers
  marker.celltype$Cell_type <- marker.celltype$gene
  marker.celltype$Cell_type <- with(reference, celltype[match(marker.celltype$Cell_type, gene)])
  return(marker.celltype)
}

Many variants of this reference table and function can be created, feel free to modify it.

In the future, I think we can expect such classifier software that uses reference data from the Human Protein Atlas, HCA and similar projects.


Feature extraction methods that can handle inconsistent numbers of atoms for molecular dynamics

This answer is no longer on topic. The question is about what is the best machine learning algorithm to use to analyse atomic coordinates from an MD trajectory in order to infer a novel properties.

  Also, it should be preface for any confused reader that unfortunately "feature" is a technical term in both genetics and statistics. For the former a feature is a region of DNA, RNA or protein that has a given property, such as open reading frame or PFAM domain. While in the latter a feature is a preferably independent property that forms part of a feature vector for different types of machine learning.

Heavy atoms are common

Assuming the pH difference only results in protonations, the protons differ between the two conditions, but not the heavy atoms. These will retain all the side chain information.

Merit of backbone only

(Section no longer relevant due to OP edit)

The common metric used to compare two structures is root mean square deviation (RMSD) of the alpha carbons to assess how much does the backbone differ and with analysis of MD simulations it is used to assess whether one is more dynamic than the other (e.g. this Bioinformatics SE question). The reason why the RMSD is done with the alpha carbons and not the minimum common subset of atoms is that backbone changes are large scale, while the side chain changes only affect locally —so their inclusion would needlessly dilute the score. The C prime and backbone nitrogens behave like the C alpha so there is no need for these.

Extraction of coordinates

Assuming by features you mean atomic coordinates, there are loads of ways. Using Pymol module for Python, one could do something like:

import pymol2
with pymol2.PyMOL() as pymol: #can be parallelised
    pymol.cmd.load('trajectory.gro', 'mytraj')
    pymol.cmd.remove('h.') #remove hydrogens
    pymol.cmd.intra_fit("name CA", 1) #align all steps in the trajectory
    for i in range(1, pymol.cmd.count_states('mytraj')):
        for atom in pymol.cmd.get_model(f'state {i}').atom:
            atom.coord #List[int int int]

One could change pymol.cmd.get_model(f'state {i}') to select state {i} and name CA for Cαs only.

Edit: Mutated residues

In the case of mutated residues or residue differences due to the two proteins being close homologues, "mutating" the differing residues to alanine, i.e. by stripping all heavy atom beyond the Cβ (keeping C′, Cα, H, N, O & Cβ) in the MD trajectory —or glycine if one is glycine. If they are homologues with different lengths, this approach would mean skipping residues that match up with gaps in the MSA. Even though not uncommon for statistically analyses, it does pose some issues statistically. It may unbury hidden residues, but length variant generally happens with loops.

In the snippet above, the pymol mutagenesis wizard could be used, but that is overkill, does not work in parallel and the approach does not port to other PDB reading modules in whatever language one wants. So say you have a list to_be_trimmed of integers of PDB residue indices (not pose residue indices):

if atom.resi not in to_be_trimmed or atom.name.strip() in ('C', 'CA', 'H', 'N', 'O', 'CB'):
      do whatever with atom.coord

In acidic environments, you get some funky things happening like exposed aspartate-glycine residue pairs forming isoaspartate or exposed aspartate-proline resulting in a cleavage, etc. For these cases, which are nearly never modelled, ignoring these values is the best option.

Edit: small comment on SASA

The investigated property list is very extensive and some reader may adopt them. So I thought I should mention that per residue SASA is generally converted to relative surface area, which is residue independent, by dividing it by the maximum SASA for that residue (values from Tien et al. 2013: 'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}), where RSA > 0.2 is an exposed residue.


Is there any value in scaffolding the output contigs of MEGAHIT assembler given a metagenomic dataset?

Update 2: It looks like your approach has actually been suggested here as one way to use the PE information. I guess megahit may not be really using the PE information anyways. I still believe that it's kind of weird but other people do suggest it, so maybe it's worth trying what Torsten suggests.

I think that scaffolding with PE information from the reads used as input to Megahit is somewhat sketchy. You could still try it, but I'd be worried about artifacts, just because de novo assemblers are so heuristic.

However, I think that it is perfectly ok to scaffold using orthogonal data. Here are some examples of orthogonal data:

A reference assembly (weird for metagenomics but still possibly ok).
long reads, e.g. ONT or PacBio (low coverage probably fine).
Proximity ligation data, e.g. Hi-C reads.
Mate-pair reads (also lead to artifacts at a somewhat high rate, but technically valid).
Linked-reads (e.g. 10X).
Probably other methods as well.

The tools used in each case would be somewhat different. For a little more information about how you might use this, here is a recent review of tech.

Full disclosure: I work for a company that sells Hi-C kits for such applications.

Update: Realized that I missed one part of the Q. I think that visually exploring the assembly using e.g. bandage is always a good idea. Quite possibly you can make some scaffolding decisions that way, but to me it sounds somewhat painful to do, especially in a metagenome where there are going to be a lot of multiple-branching collapsed regions.


how to change the UMAP use in the dimplot and feature plot

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'


Populations genetics and dynamics of bacteria on a Graph

I assume you are simulating a null distribution. Are you investigation recombination??
My main advice is to use population genetic terminology rather than geomometry to describe your simulation (e.g. vertex is a migration event between allopatric populations). Migration is investigated between populations investigated via $F_{ST}$. You appear to be assuming clonality, and population movement appear to follow a 'stepping stone model' (I might have misunderstand the question).
Mutation rate ... I'll have a guess at $10^5$ per site per year. Migration rate depends on the bacterial species being modeled, eg. MRSA spread in a very similiar fashion to that described in your model (along a moterway), in contrast soil bacteria don't adhere to allopatry at all (when I looked at them).
Fixing one parameter (population size) whilst investigation a second (migration) would be normal to assess the parameter space initially. Bacterial populations do undergo selective sweeps and using an appropriate growth model I would look at a standard population growth model (e.g. the one Bob May investigated in chaotic interactions within a differential equation). There are several population growth models and fixing migration whilst assessing the model would seem sensible. The parameter space in the combined model of migration, space and population size becomes complex.

Okay I get your model now. Its fine, the thing I'd be cautious about is the population density of the bacterial host because this would skew migration across the nodes of your graph. A quick literature review will demonstrate the importance,  double paeroto log normal distribution and its cousin (the name escapes me .... ) are sort of key words that might flag up. Bacteria are not my active thing so I'm rusty.


Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via #!/bin/sh or e.g. #!/bin/bash if it matters), so won't be affected by your choice of user shell.

If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell.  See Csh Programming Considered Harmful and other essays/polemics.

A Bourne-style shell is pretty much the industry standard, and if you choose a substantially different shell you'll have to translate some of the documentation of your bioinformatics tools.  It's not uncommon to have things like

  Set some variables pointing at reference data and add the script to your PATH to run it:

export FOO_REF=/path/to/stuff
export PATH=/path/to/foo-x.y:$PATH
foo blah blah

These will typically be shown in Bourne-shell syntax.  By using a different shell you have to translate the export commands to your local syntax, and especially PATH munging is somewhat shell-dependent.

If you're experienced in Unix, this will be only a minor niggle.  If you're a beginner, IMHO this will add a non-negligible amount of friction on top of all the other things you're learning.


How to extract RNA sequence and secondary structure restrains from a PDB file

I suggest you take a look at rna-pdb-tools we do way more than you need! :-) The tools can get you a sequence, secondary structure and much more using various algorithms, and all is well documented http://rna-pdb-tools.readthedocs.io/en/latest/

To get sequence http://rna-pdb-tools.readthedocs.io/en/latest/main.html#get-sequence

$ rna_pdb_tools.py --get_seq 5_solution_1.pdb
> 5_solution_1.pdb A:1-576
CAUCCGGUAUCCCAAGACAAUCUCGGGUUGGGUUGGGAAGUAUCAUGGCUAAUCACCAUGAUGCAAUCGGGUUGAACACUUAAUUGGGUUAAAACGGUGGGGGACGAUCCCGUAACAUCCGUCCUAACGGCGACAGACUGCACGGCCCUGCCUCAGGUGUGUCCAAUGAACAGUCGUUCCGAAAGGAAG

or you can get sequence and secondary structure via x3dna (http://x3dna.org/)

[mm] py3dna$ git:(master) ✗ ./rna_x3dna.py test_data/1xjr.pdb
test_data/1xjr.pdb
 >1xjr nts=47 [1xjr] -- secondary structure derived by DSSR
 gGAGUUCACCGAGGCCACGCGGAGUACGAUCGAGGGUACAGUGAAUU
 ..(((((((...((((.((((.....))..))..))).).)))))))

The problem is not trivial. I hope rna-pdb-tools has stuff that can solve this issue for you (http://rna-pdb-tools.readthedocs.io/en/latest/want.html#module-rna_pdb_tools.utils.rna_x3dna.rna_x3dna)


Selecting sites from VCF which have an alt AD > 10

using vcfilterjs

and the following script:

function accept(vc)
    {
    var i,j;
    for(i=0;i< vc.getNSamples();++i)
        {
        var genotype = vc.getGenotype(i);
        if(!genotype.hasAD()) continue;
        var ad = genotype.getAD();
        /* loop over AD starting from '1' =  first ALT */
        for(j=1  ;j< ad.length ;++j)
            {
            if(ad[j]>10) return true;
            }
        }
    return false;
    }
accept(variant);

usage:

java -jar dist/vcffilterjs.jar -f script.js Test.vcf


What does it mean for an amino acid to match with itself 99% of the time in a PAM1 matrix?

It means that the 99% of conservation of amino acids is just a mean value for all the amino acids. That there can be amino acids that are less conserved and some that are more conserved.

So from the 21 basic amino acids, arginine is more conserved than lysine in those alignments. But in general they are conserved the 99% of times.


Is there a difference between any scRNA-seq visualization and pseudotime?

Pseudotime is the order of cells along a trajectory which in the simplest case is a minimum-spanning tree based on a (low-dimensional) manifold such as PCA or UMAP. Without a defined trajectory you have no pseudotime, and therefore UMAP is not pseudotime, but in fact can be used to define the trajectory.


10X cellranger count, [error] The chemistry was unable to be automatically determined

The issue was that I supplied wrong indexes during fastq files generation, however, they were generated successfully, and only the next step failed with not that informative error.


How to define an outgroup to build a robust amino acid tree

I think you have a bug

raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /Users)/username/Desktop/alignment_file.phy

Try the above, where username is your username. I'll check my codes later. The most important bit is the space between the -p and word outgroup. You can omit the -p 10000 it is only needed  for very geeky maximum likelihood and sets the random number stream.


lower mapping rates in salmon v0.13 compared to previous versions

[Most likely answer based on comments attached to the question]
Hi Courtney, We've looked into the data and identified the source of the different mapping rates. Specifically, the cause is discordant reads; in this case reads that are dovetailing with respect to their positions on the reference. If you run Bowtie2 with "RNA-seq" flags (i.e. disallowing discordant mappings), then you get essentially the same mapping rate as salmon here. Is allowing dovetailing reads something that should be very important in your analysis? If so, we can consider what would be necessary to add back support, but generally, they provide evidence of tenuous quality.


Can foldx position-scan produce different result with different lists of mutations?

Minute values.
So just a reminder that the ∆∆G units is a kcal/mol. Water bounces around and collides with different energies following a Maxwell-Boltzmann distribution, where the average collision at 25°C is 0.59 kcal/mol (thermal energy constant).
A hydrogen bond reduced the ∆G by 1–3 kcal/mol. Other interactions such as π-π, π-S are on the 1–2 mark. So the difference you are seeing is totally negligible.

Heuristic.
The reason why you see the differences is that the force field calculations are heuristic and are not testing all perturbations —it is not only sampling a limited amount of conformers of side chains (kind of like PyMOL would be if it were FF based), but neighbour perturbations. So sometimes you hit a wee energy minima. Hence why in some protocols replicates are required —not mutation scanning mind you.

Bonus: Protein pI.
Lastly, a caveat about pI of protein, you probably already know this, but best be safe. Your scores are all low, so I am guessing these surface mutations? Improving thermal stability is great for increasing the T_M measured by DSC or catalysis if you are scanning a transition state as done for theozyme design, but it does not mean you get better solubility. In this case, you are introducing positive charges, which will alter the pI of the protein upwards and most protein have a pI less than 7 and the closer it is to pH 7 the more they crash out. If you are designing a high pI GFP, Rosetta supercharge does FF calculations to alter the pI of protein, which might be worth a gander.


Is there a computational tool or possibility to identify mRNA isoforms from the count matrix of a bulk RNA sequencing dataset?

You can back to your bam file and use cuffdiff tools it is very able to do your job follow this to liks
http://cole-trapnell-lab.github.io/cufflinks/manual/

https://github.com/Jeanielmj/bioinformatics-workshop/wiki/Differential-Analysis-with-Cuffdiff

My best


Mates in a pair are giving signal for a single strand, though this is only apparent in bisulfite-treated samples. The reason for this is that you are sequencing the ends of fragments of material and it's random which strand of this (or both in different clusters) will get sequenced. Different pairs of reads will then randomly sequence one strand or the other.

Single-end reads are the same as paired-end reads, the second read in the pair is just not made.