Merge multiple PDB files into a single file I have several pdb files and each of them corresponds to a domain on a protein. I want to merge these files into one single pdb for the whole protein, and the approach I'm looking for should allow me to order these domains according to the order by which they're located on the protein (for example, merging domain A, B, and C by the order CBA). It would be preferable if I can do the merging with VMD but any suggestion would be appreciated.

PyMOL is a rather flexible tool for doing this.

The command create combo, domA or domB or domC in the PyMOL app will create a single object. However, they will need to be the same chain and segment (segi). Also you will need to fix the numbering to match what you expect them to be.
To fix this you need to do:

 alter domA, chain='A'
 alter domA, resv+=0
 alter domB, chain='A'
 alter domA, resv+=100
 alter domC, chain='A'
 alter domA, resv+=200
 alter all, segi=''
 sort
 create combo, domA or domB or domC
 save combo.pdb, combo

PyMOL can be installed as a Python 3 module, which works in the same way but with pymol.cmd and is a lot more replicable. But for starting out the App is fine.

import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.alter('domA','chain="A"')
     etc.

Now. This assumed that the coordinates are where we want them to be. Which is highly unlikely and that there are zero overlaps.

If there are overlaps, aligning the overlapping residues (align domB and resi 80-100, domA and resi 80-100) can be done. Just remember that the first element is the mobile object, the second the target, hence why I wrote it backwards.
If do not have overlaps and want the domains docked together you can use the DEMO tool by the I-Tasser people.
If you want your domains as "balls on a string" aligned on the x-axis for a nice picture and you do not care if it may not be biological, you may want to check out a Github repo of mine for doing just so (it offsets, aligns if possible and spaces them out, but does not add linkers).

If you want to add linkers, you can use submit your sequence to Phyre or I-Tasser with your missing density model as a reference, or you can add them with rosetta remodel (tricky) or make them up in Coot or PyMOL (bad idea, but if you are in a hurry... here is a tutorial of mine).


What are the units for the force output by nstfout, the velocity output by nstvout, and the coordinates output by nstxout in gromacs? What are the units of force, velocity, and coordinates output by gromacs trajectory? In .mdp, they are specified by nstfout,nstvout,nstxout.

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


OK. This is a bit convoluted, but it works. There might be simpler solutions...

First we'll write a function that takes a data frame and a character position (ie; first or second position in the string), groups the data frame by the second and third columns, counts occurrences of each character in the specified position of the 18th column, and returns a vector according to your ranking scheme (this doesn't include counts of N since it doesn't occur in your example data, but it's a simple extension of the case_when line to add it):

library(dplyr)
library(tidyr)
library(stringr)

summarise_characters <- function(df, pos) {
  df %>% mutate(character=str_sub(X18, pos, pos)) %>%
  group_by(X2, X3, character) %>%
  summarise(n=n()) %>%
  ungroup() %>%
  spread(character, n, fill=0) %>%
  mutate(pos1=case_when(`@` > ß ~ '@',
                        `@` == ß ~ 'N',
                        `@` < ß ~ 'ß')) %>%
  `$`(pos1)
}

Next we remove duplicate lines from the data frame, then create a new column, pasting together the output of our function called on the first and second character:

df %>% group_by(X2, X3) %>%
  slice(1) %>%
  ungroup() %>%
  mutate(counts=paste0(summarise_characters(df, 1), summarise_characters(df, 2))) %>%
  select(X1, X2, X3, X4, X11, X12, X13, counts)

# A tibble: 3 x 8
     X1        X2        X3                   X4       X11                  X12       X13 counts
  <int>     <int>     <int>                <chr>     <int>                <chr>     <int>  <chr>
1     1  90042011  90042031 AGTGCCACCAGGTGGGGCCG  90184354 GAGGCCACAAGAGGGCACAA  90184337     @N
2     1  94853396  94853416 GCCGCCACTAGATGGTGCTA  94969274 CAGTCCTCTAGAGGGAGCCC  94969266     ß@
3     2 103584283 103584303 TGGGCCACTAGGAGGCACTG 103841456 CATGCCACAAGAGGGCATCA 103841456     NN


What's the scaling for HOMER metagenes? I'm trying to use HOMER to make a metagene profile over gene bodies using a bedgraph file I've generated. The problem is that every time I do, I get really weird scaling on the y-axis. I should be getting average values across the gene body on the order of 5-10, but instead I'm getting values on the order of 0.03-0.05 or less. The weird thing is, when I don't use metagene script or when I don't use histogram with -size given, I get perfectly normal values -- but that's not what I want, unfortunately, for the metagene profile.

What am I missing?

Figured it out! For anyone who finds this thread in the future and is wondering what is going on: HOMER, when using the -size given option, normalizes the y-axis to the number of basepairs in the intervals -- in this case, the length of the gene body. To get the "raw" values, as I wanted, you need to multiply by the length of the gene body to "undo" the normalization.


How to force pilon to correct a reference with low coverage mapped reads? I have mapped some reads to a reference and want to change SNPs in the reference to SNPs found in the reads, even at low coverage. For example, here is a screenshot of one reference in IGV that I would like to correct.

When I use pilon to try to correct the reference, however, no bases are ever changed, even if using the option --mindepth 1.

Here is the call that I am using:

java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --frags {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1

How can I use even low-coverage information to correct my reference?

Pilon will not consider fragments that are not proper pairs if the --frags option is used to pass the alignment .bam file. The reads mapped in the above IGV screenshot are unpaired.

Instead, if you want Pilon to consider all reads that map as legitimate and usable for correction, pass the .bam file with the --unpaired option instead.

To above corrected example would be:

java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --unpaired {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1


How to run Cox proportional Hazard model for each group in R? My population is divided into 10 groups. I want to run the cox-proportional hazard model separately for each of these 10 groups so that in the end i can compare which groups has highest hazard ratio depending uopn the covariates.
How can i do this in R? Earlier i was using this formula.

Model_1 <- coxph(serve_object ~ GROUP + age+ Smoking+ Diabetes, data = DF
summary(Model_1)

But the output is very confusing. it gives hazard ratio of 9 out of 10 groups. I think it dealt group as covariate.

The fact that you supplied GROUP as a covariate, you are getting 9 coefficients as one of the ten levels is used as the reference level.

If you would like to fit one model per group (and in this scenario you would not be correcting for GROUP), you should first split your data frame and fit your models on the resulting list's elements:

lapply(split(DF, DF$GROUP),
       function(x) coxph(serve_object ~ age + Smoking + Diabetes, data = DF))

Example with the mtcars data set when grouped by number of cylinders, cyl having values of 4, 6 or 8:

> lapply(split(mtcars, mtcars$cyl),
+        function(x) lm(disp ~ hp, data = x))
$`4`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
    59.0369       0.5579

$`6`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   290.9178      -0.8799

$`8`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   320.2083       0.1572


DoHeatmap Hierarchical Clustering Seurat I have a list of genes that I'd like to visualize using the DoHeatmap function in Seurat. However, the output of the heatmap does not result in hierarchical clustering and therefore makes it very difficult to interpret. My code is,

 DoHeatmap(object = obj, genes.use = genes),
          slim.col.label = TRUE, remove.key = T,group.label.rot = F, use.scaled = T)

Is there a way to adjust the DoHeatmap command to rank the cells by the intensity of gene expression? Here's an example output:

You need to order the marker matrix (e.g. by avg_logFC) before calling DoHeatMap.

library(dplyr)
all.markers <- FindAllMarkers(object = obj)
top20 <- all.markers %>% group_by(cluster) %>% top_n(20, avg_logFC)
DoHeatmap(object = obj, genes.use = top20$gene, slim.col.label = TRUE, remove.key = TRUE)


Extract sequences from partial Header I have a library I would like to split them based on their headers. Some are related to RNA and DNA.
The header contains much information, but the most important the presence of DNA or RNA/LTR.., these partial words could be in between or at the beginning.
The point is the knowledge of how to extract sequences from a partial header occur in between the ID
My question is how to use grep or awk to grep the header that has one of these words along with sequences?
Note that The sequences are more than one line.
Or perhaps grep the specific word and ignore what before and after?
>Tigger16a#DNA/TcMar-Tigger DF0000028 TcMar-Tigger **DNA** transposon
>rnd-4_family-38#SINE/MIR ( Recon Family Size = 20, Final Multiple Alignmen
>rnd-6_family-31751#LTR/Gypsy ( Recon Family Size = 26, Final Multiple Alignment Size = 22 )
>RNA2558#LTR/ERVL
>NonDNA1#LINE/I-Jockey
>DNA5#DNA/TcMar-Tc1

I have tried to create a list of the Required IDs using grep creating the list and extracting the sequences but for some reasons, the output has more sequences than those specified in the DNAID.txt list.
grep -A1000 -w -f DNAID.txt.fa MyLibrary > DNA_Sequence.fa

You could use Biopython. If its a fasta file, it will be complicated to write fasta output (especially multiline fasta) with grep or awk. Simple solution is to use biopython, so that you can even match for any complex patterns in the fasta header.
from Bio import SeqIO

rna_records  = []
dna_records = []

for seq in SeqIO.parse("in.fa","fasta"):

    if "RNA" in seq.id:
        rna_records.append(seq)

    elif "DNA" in seq.id:
        dna_records.append(seq)

SeqIO.write(rna_records, "RNA_out.fa","fasta")
SeqIO.write(dna_records, "DNA_out.fa","fasta")

This code needs to be optimised if you are dealing with large sequences like human genome fasta as it holds everything in memory.


Extract autodocked protein-ligand connections programatically I have 2 crystal structures, one for my protein and one for my ligand (I have several protein-ligand pairs). I am using AutoDock Vina to simulate docking, which returns another file of the ligand with updated 3D coordinates. Everything is in pdbqt format, but I can use OpenBabel to convert them if needed.

I am looking for a programmatic way (mainly python, but other languages work as well) to load both the protein and the docked ligand, then extract only the acting atoms which connect the 2 structures together (and the bonds between them).

If a direct answer is too complex, I would be happy if someone could suggest a tool/framework I could use to manipulate this sort of data easily so I can implement a solution on my own.

With pymol as a python module (conda install -c schrodinger pymol, not the GUI one), it is very easy.

import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load('protein.pdbqt', 'protein') # the second is the name within pymol instance.
     pymol.cmd.load('ligand.sdf', 'ligand')
     pymol.cmd.save('neighbours.pdb','ligand expand 3') # get atoms 3 Å around

Done!
However, if you are doing stuff in python in the first place you might want the last line as a string: pymol.cmd.get_pdbstr('ligand expand 3'). Changing the selection 'ligand expand 3' to 'byres ligand expand 3' will save the full residues.

Alternatively you can get the atom coordinates themselves etc. with it. To see what attributes Atom has you have you can do this:

for atom in pymol.cmd.get_model('ligand expand 3').atom:
    print(dir(atom))
    break

Which will give a large amount of attributes, such as atom.coord, which is a 3-element list.


GOseq analysis with evidence code filter When using GOseq analysis on RNA-seq data, I often find many 'false positives'. What I mean with that is that some genes in certain categories are not really involved in the process, but are only involved due to 'Electronic Annotation' (IEA). Most of the genes are IEA. A list of all evidence codes is here.

Is it possible to (easily) do GOseq analysis without this evidence group?

I mean, I know how to do it the hard way (make a custom GO-genes annotation in e.g., biomaRt removing all IEA evidence code genes, and perform GOseq with custom genome). But is there also a easy way in GOseq available?

A bit weird to answer my own question, but for who's interested I can show the 'hard' way using biomaRt and then use goseq with gene2cat defined.

So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.

library(biomaRt)

table <- read.table(file, header = T, stringsAsFactors = F)

ensembl <- useMart("ENSEMBL_MART_ENSEMBL")

ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)

allGenes <- table$Geneid

GOs <- getBM(attributes=c('ensembl_gene_id', 'go_id', 'name_1006', 'go_linkage_type', 'namespace_1003'),
      filters = 'ensembl_gene_id',
      values = allGenes,
      mart = ensembl)

dim(GOs)
[1] 477082      5

table(GOs$go_linkage_type)
          EXP    IBA     IC    IDA    IEA    IEP    IGI    IMP    IPI    ISA    ISM    ISO    ISS    NAS     ND    RCA    TAS
 36931     36  27918    525  34038 210114    630   6891  30219   6205   4234     20  77365  24673    609  12821     99   3754

# I only want BP and no IEA evidence code
GOs_no_IEA <- GOs[GOs$go_linkage_type != "IEA" & GOs$namespace_1003 == "biological_process",]

dim(GOs_no_IEA)
[1] 111944      5

Then you can continue with goseq.

library(goseq)

All <- read.table(input_All, sep = "\t", header = T)

All_genes <- as.character(All$Geneid)

TopTable <- read.table(input, sep = "\t", header = T)

sig_genes <- as.character(TopTable$Ensembl)

# Make gene vector

gene.vector <- as.integer(All_genes %in% sig_genes)

names(gene.vector) <- All_genes

gene.lengths <- All$Length

pwf <- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
head(pwf)

GO.wall <- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")


Scaffolding a genome with hybrid data I am assembling a ~500MB genome, and have ~150x long reads and ~200x 150bp PE short reads, with a ~400bp insert size.

I've done a lot of work with minimap+miniasm, and have what I think is a good set of unitigs and a .gfa file from miniasm. Depending on the parameters, I have a few hundred to ~1K unitigs. I can easily polish these unitigs with pilon, racon, and nanopolish, so that they are as error free as possible.

The question I have is: given access to

Polished (or not) unitigs
150x long reads (long enough to span almost all repeats in the genome)
200x 150bp PE short reads (~400bp insert size)

What would you do to scaffold a genome? And what are the pros and cons of your suggested approach?

What I need here is a draft genome. I would prefer a conservative approach that is more fragmented and has fewer misassemblies. So of course one option is just to keep the unitigs and make no attempt at scaffolding. But it does strike me there may be suitably conservative approaches I haven't considered. E.g. perhaps I could use ABYSS to scaffold my polished unitigs using my short-read data, thereby leveraging whatever information I can from the PE short reads?

Try LR_Gapcloser.

I've used L_RNA_Scaffolder for trying to scaffold a genome (which turned out to be more complex than I had expected). It looks like LR_Gapcloser (written by the same people) is similar, but designed specifically for scaffolding using long-reads.

That page also suggests PBJelly and GMCloser as competing tools.


Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type? Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type? For example, a FASTA file has 2 sequences. Can the first one encode amino acids while the second one encodes bases?

Thank you

While there's nothing stopping anyone from doing that with the FASTA format (after all, it's just a text file with '>' defining header lines), I don't know of any software that would support such a file structure. At best, it would interpret the nucleotide sequences as protein sequences (A/C/G/T are all valid 1-letter protein codes).

A better question to ask would be "Does ultra-cool-bioinformatics-tool X support combined nucleotide and protein sequences in the same FASTA file?" In which case the answer would most likely be, "No."


How can the Autocovariances, autocorrelations, and autocorrelation coefficients be calculated from a Protein Amino Acid Sequence? Given a normal protein sequence with the 20 standard amino acids, how can the 'Autocovariances', 'autocorrelations', and 'autocorrelation coefficients' of the sequence be calculated?

What is meant by these terms in the context of a protein sequence?

What are the appropriate algorithms and formulas?

Is there any existing package/function to calculate these values?

Protein autocorrelation is similar to usual autocorrelation of time series (that is, function $r(k)$ measuring correlation between values lagged by $k$), except that instead of time dimension we slide along the backbone of amino acid chain. The "dependent variable" can be any quantifiable property of the amino acids - then the autocorrelation measures similarity in that particular property.

See e.g. Fernandez et al., 2008, for more extensive background:
https://www.tandfonline.com/doi/full/10.1080/08927020802301920?scroll=top&needAccess=true


Finding orthologues using BLAST on the NCBI database I'm an informatics student who has essentially zero knowledge of biology. I BLASTed my gene and have 1000s of results with very low E values. Where do I go from here if I want to find orthologues?

BLAST is often a sensible way to find gene/protein homologs, but resolving those as orthologs vs paralogs etc. is a non-trivial task, especially when considering large numbers of species.

Since you are new to bioinformatics I recommend you search an ortholog database that already exists.

Here are some databases:

orthoDB
OMA
EggNog

All of these should allow you to search either with text or a sequence to identify the orthologous groups that your input belongs to at different taxonomic levels. Sequence searches either follow pairwise alignment (e.g. BLAST, DIAMOND, RAPSEARCH) or profile HMM type (e.g. HMMER3) approaches.

A more complete list of ortholog databases is available from "Quest for orthologs".


Fastqc- Per Base Sequence Quality I am trying to figure out how to interpret the "Per Base Sequence Quality"? What does Position in reads (bp) mean? Also in order to draw this box-plot graphics, more than one quality scores are needed so we can sort them and calculate the mean, median and quartiles. In this diagram, how these values were calculated for each position?

There is a pretty good description here. The position in the reads means literally the position in the reads, so position 2 means the quality of all bases on position 2 in the fastq file (all reads). In the description link, you can find how the boxplots and blue and red lines are calculated (and what they mean). They are all based on the total quality values of all bases on a certain position, all bases means all bases of all the reads on a particular position.


Ensembl protein identifiers from different assemblies In order to obtain orthologs for human genes, I am working with InParanoid version 7.0. This version of the database uses Ensembl Protein Identifiers (ENSP) from ENSEMBL version v54, which is based on NCBI 36 assembly of the human genome.

For my project, I need to have the Ensembl Gene Identifiers (ENSG) of these proteins, but in the last version of ENSEMBL, which is based on GRCh38 assembly. I need to use this version of InParanoid for consistency with previous analyzes.

Thus, my question is, could I just use biomaRt to map ENSP (from NCBI36) to ENSG (actual version) or this approach is inherently wrong?

Is it possible to "convert" ENSP identifiers from NCBI36 to GRCh38 and then do the map to ENSG?. I am aware of liftOver but I am not sure it would work here.

Those IDs are elderly! Ensembl 54 was 2009!

I would recommend using BioMart combined with the ID history converter. The ID history converter will convert old IDs to new, and BioMart will convert ENSPs to ENSGs. You can either use the ID history converter with the ENSPs first, then the current BioMart to get the ENSGs. Or you can use Ensembl 54 BioMart to convert the old ENSPs to old ENSGs, then update the ENSGs with the ID history converter. I would probably do both and compare results.


differential gene expression complex design no replicates We have an experimental design as seen below

Where we administered drug at 0 min for each mouse genotype and took them down at given below intervals. wt and ko mouse models  were administered only with the vehicle which would be our negative control.

Performed rnaseq (total rna) truseq stranded protocol would like to see differences between each time point compared to the vehicle, used star aligner to align to mm9 , htseq for gene counts .

sample  type    time
WT  Vehicle     30 min
KO  Vehicle     30 min
WT  Drug    30min
KO  Drug    30min
WT  Drug    1 hour
KO  Drug    1 hour
WT  Drug    2 hour
KO  Drug    2hour
WT  Drug    3 hour
KO  Drug    3hour

We do not have any replicates and hence cannot use any of the differential expression statistical tools (asked this question at bioconductor where Dr. Love himself replied link to question asked) where he suggested a linear model.

we were wondering if someone could help us with any such design (papers/previously performed analysis) what we would like to see is comparison/differences (significant) between these groups
(veh vs 30min, veh vs 1 hour, beh vs 2 hours, veh vs 3 hours)

Although it is not recommended to use no replicates, in the edgeR manual they give some advice on how to go on with no replicates design. See page 21 of user guide. You can e.g. estimate a BCV value. Of course it is still a trick and not sound statistics.


Problems to Extract uniquely mapping reads from BWA MEM alignment I did a mapping of genomic paired-end reads to a reference assembly using bwa mem. I need to extract the reads that mapped only once to my reference.
For that, I have tried to follow the method proposed in another post Obtaining uniquely mapped reads from BWA mem alignment and unfortunately in my case, after using:

$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d

I still obtain some outputs. Here is an example:

NB501461:7:HH3GGAFXX:4:21612:26866:20400    99  ctg7180000256968    441 0   147M4S  =   560 230 AAGCTTTAGTCGCATGGCCATTCACAAGACCAGTGATTTTTTATAGAAGTCGTTCAAGAGATACTGATGAGGAGTACAGCCATTCTTTAAGCTATACAAAAAGCACCACCCACACCACACAATCCTCAAAAAGATCAAGGCAGTACAGATT AAAAAEEAEEAEAAEEEAAEEEAEEEEE/EEEEEEAA/EEEEEEEEEEEEEEAEE/EEEE///E6EAEEEEAEA<EE/EA/</EEAE/EEEEE//EEEEEE6E/<EAE/EEE/EE</E/E//AE<6EA//A<A///<A/6A/6//A///// NM:i:7  MD:Z:34A56T11A4A3A2A22T8    AS:i:112    XS:i:117
NB501461:7:HH3GGAFXX:4:21612:26866:20400    147 ctg7180000256968    560 0   32S111M8S   =   441 -230    CCAGCTCCACCACCCGACCCCCCCACCCACCCCAATCCTCAAACAGATCCATGCAGTACATTTTCATTTGTTCTGGCCCTGCGGCTCTTAGTTTAGGAAATATTCAATGACTTGCTACCTAGTAAGAACATTATTGTGCACAACGAATAGA ///E/////A////A/6A/EE/<//////A//AA/EEE/A/<E///A/A<//</////A/////<///EE////E/E/E/E</E//EA//E//A/</E//E//AE<//A<///EE///////E//E//EA/A//EA/AE/A/E/E//6/A/ NM:i:7  MD:Z:11A5A22A4A4T27A7A24    AS:i:76 XS:i:79
NB501461:7:HH3GGAFXX:4:21612:6399:20400 113 ctg7180000282074    357 8   73S50M28S   ctg7180000432624    28  0   CTTATCCAATTTTTAAACATAGGATTAGGGTACCTTCCAAATGTTATTGAAACTTAATCTGAAGATGACCTCAAATGGAAGATGACCTCCAATGGAAGACAACCTTCCATGGAAGACGACCTCAAAAATGAAGAAGAAGAAAACATCTAAG /A//<//<<//E/66A//<A////A/<</E/EA/EE//<E<E//<AE//A/EEAEEEE/E//A/A/</EEE//E<EAEEEEE/EEEEEE<EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEAAAAA NM:i:2  MD:Z:32C1A15    AS:i:40 XS:i:35
NB501461:7:HH3GGAFXX:4:21612:6399:20400 177 ctg7180000432624    28  26  76S75M  ctg7180000282074    357 0   CTCCTTACTGTTGTTTCTTCTTTTTTCTTCTCTTTCATTTTTGAGGTTGTCTTCCATGTTTTGTTGTCTTCCATTTGAGGTCGTCTTCATATTAAGTCTCTATAACATTTGGCTTTGTCCCCTAAGCCTAAGTTTCAAATTTGGAAAAGAT /////</////6/////A////EAEE<//////<A////E/E<//A//<///E//<E/////E/A/E//E/EAEEAE/EA///E//E/AE/<EAEE/6E6EEA/EE/AEEAE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA NM:i:4  MD:Z:5T15T17G29G5   AS:i:55 XS:i:44

As I do not see any XA or SA flags there anymore, I wonder what those can be?

I have tried to continue the steps I want to pursue until the "ValidateSamFile" from Picard tools, and there I do not know if it is linked or not, but I get the message:
"Error: mate not found" for a few thousands of reads.
Any idea what could be wrong, and how to fix it?

I hope I was clear enough, otherwise please feel free to ask me for specific details.
Any help is welcome,
Thanks!

The issue is that you are filtering reads without regard to the pairing so you might filter out one of the mates if one aligns uniquely and the other does not.  You will need to filter out the un-mated reads from the file to have a valid BAM.

  As I do not see any XA or SA flags there anymore, I wonder what those can be?

I am not sure what this means.  The grep should remove those flags and only leave records without them.

As a side note if you don't want split alignments I would recommend using a different aligner such as BWA aln or BowTie2


TEs vertical or horizontal transfers? After predicting and annotating all the transposable elements. For four related species.

Now I would like to find out these predicted candidates found in my costume library, either vertical or horizontal transfers?

are there any methodologies to suggest?

You make a phylogenetic tree of the "host species" and compare it to the phylogeny of the TE's. Homology is not a good idea for evolutionary trees.

There are two snags:

Concerted evolution, the TE's will risk recombining with each other within a host genome and that can mess up the phylogeneic signal (the theory as why is quite complicated)
If host species phylogeny are higher eukaryotes thats a lot of genetic material, so you'd need a short-cut. You can either use someones species tree, or focus on TE's within a subsection of the host genome and produce a tree of this subsection.

With point number 2 "a subsection of the genome" you can then assess synteny, i.e. has the TE jumped between position between species. If the TE hasn't moved then its an orthogonal proof of putative vertical transmission.

What you are looking for is incongruence between host phylogeny and TE phylogeny as evidence of horizontal transmission and congruence to call vertical transmission. However, incongruence could occur due to the complexity of concerted evolution, so synteny would assist again as an orthogonal proof. Congruence is assessed as a "mirror tree", there is quite a lot of theory on this looking at host-parasite coevolution

For making trees RAxML or PhyML are flavours of the month particularly for TEs.

A state of the art-of-the-art review is here, but it goes into significant detail into precise proofs, which in my opinion arn't needed


how to sort and join two files based on First Column Id I have got two files large files like this in Tab-delimited format, trying to merge in R, do we need to sort the file before merging?
 otu <- read.table("asv.txt", sep="\t", header=T, blank.lines.skip=F, check.names=F)
> head(otu)
                ID                   Soil.9.S42 Soil-27.S33.L001 Soil-45.S54.L001
1 d2ec9f3b77975c0f457e4b7413b217ff         34               51               84
2 3147790f0d5a78316fb9dd64f53b9473         28               32               52
3 97aecc1f35cc1f50db507ad71dd22367         33                0                7
4 bfad6370d28182cc6304844e9bec7fb6        545              237              127
5 5fa2a987221a1d9ca416148570c18086          0                0                0
6 8e60d301122d7aa359eb6b0b00f37f62        124              143              128

> taxa <- read.table("taxonomy.txt", sep="\t", header=F ,stringsAsFactors=F,quote="")
> head(taxa)
                  ID                     V2             V3                  V4
d2ec9f3b77975c0f457e4b7413b217ff Bacteria  Cyanobacteria    Oxyphotobacteria
3147790f0d5a78316fb9dd64f53b9473 Bacteria Proteobacteria Alphaproteobacteria
97aecc1f35cc1f50db507ad71dd22367 Bacteria Proteobacteria Gammaproteobacteria
bfad6370d28182cc6304844e9bec7fb6 Bacteria Proteobacteria Alphaproteobacteria
5fa2a987221a1d9ca416148570c18086 Bacteria Proteobacteria Alphaproteobacteria
8e60d301122d7aa359eb6b0b00f37f62 Bacteria Proteobacteria Deltaproteobacteria

Now having issue with sorting file but merge command work perfectly;
newotu <- otu[order(ID),]
Error in order(ID) : object 'ID' not found
> total <- merge(otu,taxa,by="ID")

Thanks for your help.

You should check left_join() method in dplyr package. In this way you do not need to sort the data frames, as the function automatically take care of this.
For example: merged_df <- left_join(otu, taxa, by = "ID").


How to find possible viral subsequences in a Eukaryote genome? I have an assembly of an algae and want to find if it has anything analogous to endogenous retroviruses' subsequences in its DNA.
My guess is I should manually blastn all subsequences until I find (or not something) that could have come from a virus.
Is there a better way to do this?

Your guess sounds reasonable to me but you should consider using KRAKEN2. You can create your custom database (you would have done it with makeblastdb anyway) and kraken will check k-mers against the database to find homologous sequences.


"megahit: command not found" after correct installation with miniconda I want to use MEGAHIT (git: https://github.com/voutcn/megahit) to assemble metagenome shotgun sequencing reads (illumina).

I install the package with the command conda install -c bioconda megahit as described in on the github page which yields:

Solving environment: done

## Package Plan ##

  environment location: /home/datura-workstation/miniconda3

  added / updated specs:
    - megahit

The following NEW packages will be INSTALLED:

  megahit            bioconda/linux-64::megahit-1.2.9-h8b12597_0

Proceed ([y]/n)?

Where I click yes and the packages is installed. However, after installation I am not able to use the command megahit as this results in the error: megahit: command not found.

I use Ubuntu 18.04.3 LTS and have installed miniconda3.

I reproduced the installation of megahit on a windows computer with a virtual machine using the same Ubuntu version. On my virtual Ubuntu I CAN use the megahit command.

Could anyone please help me solve this problem? I searched the internet but cannot find a good solution on the issue (also not on the github issue page of this package).

type conda info --envs.  Look at your list of environments.  Say it's base that has megahit in it, do conda activate base, then just do megahit.


R: 'Matrix' can not be unloaded, but 'writeMM' method not found I want to save a sparse and very large dgcMatrix onto disc. I read that it could be done with writeMM method. So, when I am trying to do that:

writeMM(UMI_count, "/home/gene_count_filtered/filtered_dataset")

Rstudio is giving me an error:

Error in writeMM(UMI_count, "/home/gene_count_filtered/filtered_dataset") :
could not find function "writeMM"

Ok, I read that this function is in the Matrix package. I installed it, but when I try loading it:

library(Matrix)

I am getting another error due to other libraries using Matrix already on the lower level, as I understand:

Error in value[[3L]](cond) :
Package ‘Matrix’ version 1.2.11 cannot be unloaded:
Error in unloadNamespace(package) : namespace ‘Matrix’ is imported by
‘scran’, ‘scater’, ‘SummarizedExperiment’ so cannot be unloaded

I also tried justwrite.matrix(UMI_count, "/home/nikita/filtered_dataset"), but got an error:

  Cholmod error 'out of memory' at file ../Core/cholmod_memory.c, line 147

What could be the solution to the problem? I just need to somehow save a very large dgcMatrix on disc.

Ultimately I just restarted Rstudio as Konrad Rudolph suggested, imported Matrix library, and then writeMM method worked just fine. However, rownames were converted into just numbers, but I needed to keep the gene names instead as it was in dgcMatrix, so after some search I found a different way of doing that from here:

http://www.sthda.com/english/wiki/saving-data-into-r-data-format-rds-and-rdata

So, I just did save() (and then, to load back - load()) and the rownames are maintained like I wanted it to be. So, storing it now as .RData file.


Importing GFF file with Biopython Is there a way to import a GFF file for an organism with Biopython in the same way you can for a Genbank file?

For example,

from Bio import Entrez as ez
ez.email = '...'
handle = ez.efetch(db='gene', id='AE015451.2', rettype='genbank', retmode='xml')
h = handle.read()
print(h)

This will import a Gebank xml file for Pseudomonas putida.  Is there a way to do this for a GFF file?

No, there is currently no GFF support in biopython.

However, you can read in GFF files into python using this package, gffutils. There are also a few other packages to read/write GFF files, like gff3.


How to demultiplex a mix of single-indexed and dual-indexed samples The problem

If I have a sample sheet that contains both single-indexed and dual-indexed samples, I can split it up into two sample sheets and then run bcl2fastq on each one. However, when doing this, large Undetermined fastq files are generated. E.g., when processing the single-indexed samples, all the dual-indexed samples go to Undetermined. And when processing the dual-indexed samples, the single-indexed samples go to Undetermined.

Additionally, because they are being processed separately, if any single index A is part of a dual index A+B, then when processing the single indexes, an A+B may be mistaken for an A, so it would seem that they need to processed simultaneously to avoid this mis-assignment.

The question

Given a sample sheet and directories of BCL files, how can such a set of sequencing data be demultiplexed correctly, either using bcl2fastq or the Picard tools?

To put it another way, I want to demultiplex a single sequencing run that contains both single-indexed and dual-indexed samples. It can be assumed that the indexes are sufficiently distinct such that any sample's index configuration is different from any other. But assuming that the different index configurations are not segregated to particular lanes of the sequencer, the question is how to demultiplex the files correctly such that both the single-indexed and dual-indexed samples are recognized.

Attempts at a solution

Using bcl2fastq directly

If a sample sheet of a mix of indexes is given to bcl2fastq (v2.17) directly, it produces the error

ERROR: bcl2fastq::common::Exception: Success (0): .../bcl2fastq2/src/cxx/lib/layout/BarcodeCollisionDetector.cpp(127):
Throw in function void bcl2fastq::layout::BarcodeCollisionDetector::validateNewBarcodeSizesAgainstExisting(const std::vector<long unsigned int>&) const
Dynamic exception type: boost::exception_detail::clone_impl<bcl2fastq::layout::BarcodeCollisionError>
std::exception::what: Barcodes have an unequal number of components.

  Barcodes have an unequal number of components.

Using picard

It seems that it should be possible using the Picard tools, but I have not found a way to set up the inputs to ExtractIlluminaBarcodes and IlluminaBasecallsToFastq that processes this configuration correctly.

The syntax for supporting multiple index configurations is not entirely clear. But when using various combinations of N and '*' on multiplex_params.tsv and barcodes.txt required by the Picard tools, a large Undetermined file is still produced and actual sample fastq files are tiny, indicating it is not processing them correctly.

Using multiple calls to bcl2fastq

As indicated in the discussion above, this suffers from the problem of having all the dual indexed samples going to the "Undetermined" file when processing the single-indexed samples, or vice versa, and creates two output directories of reports and stats, which must be merged.

Padding empty indexes with N

By padding, I mean converting ACGTACGT to ACGTACGT+NNNNNNNN so that the single-index samples in the sample sheet are "dual" as well. This strategy seems that it would work except there is a bug in bcl2fastq that it treats "N" literally instead of as a wildcard. See the release notes for details.

I'm late to the party but simply run cellranger mkfastq twice, with different arguments for the mask and with --filter-dual-index for the double-indexed samples. In a second step you then have to disentangle the dual A+B reads from single-index A reads based on the readID's that are found in both of the demultiplexings. This is cumbersome and time+diskspace consuming, but it can be done. I would use some Unix tools for that but there are probably smarter ways.
EDIT: I mistakenly said bcl2fastq, this must be cellranger mkfastq. (thanks AmadeusDrZaius)


Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter? Disclaimer: This is a self answered question for documentation purpose and I adapted this from the following github gist. Especially from users terrycojones and peterjc as well as sujaikumar who raised the issue.

I have a strange situation. I have two BLAST databases (20170821 and 20170824) that both contain a certain subject sequence LC074724. When I query them both with a given query, Q, in one case (20170824) there is a very strong match and in the other (20170821) no match at all.

Why is there not even a lower quality match found although the databases are nearly the same?

BLAST version

$ blastn -version
blastn: 2.6.0+
 Package: blast 2.6.0, build May 11 2017 22:22:40

Both databases contain an identical subject sequence LC074724:

$ blastdbcmd -entry LC074724 -db 20170824 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

$ blastdbcmd -entry LC074724 -db 20170821 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

Only difference for the databases is 20170821 is just 20170824 with an extra 643 sequences. Rest is completely identical.

I'm trying to figure out why. This will be a longer read, Tl;dr at the end:

Doing a match

There is a good match of Q against LC074724 in 20170824:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
Q   LC074724.1  2581

Note that the bitscore of the LC074724 match is 2581.

But there is no match of LC074724 at all in 20170821:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724

The option -outfmt '6 qaccver saccver bitscore' tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.

What are the best hits for 20170821

So what are the bitscores of the Q matches in 20170821?

$ blastn -outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head
3301
2616
2590
2590
2563
2554
2547
2544
2540
2540

The bitscore of Q against LC074724 in 20170824 (2581) would have placed it in 5th position in the top 10 matches of Q, but LC074724 is not matched at all in 20170821! The lowest bitscores matched in 20170821 is 2354.

What's the difference between the databases?

Get the ids of the subjects in each database

$ blastdbcmd -entry all -db 20170824 | egrep '^>' | cut -c2- | sort > 20170824.ids
$ blastdbcmd -entry all -db 20170821 | egrep '^>' | cut -c2- | sort > 20170821.ids

$ wc -l 2017082?.ids
    9754 20170821.ids
    9111 20170824.ids
   18865 total

20170821 has 643 sequences that are not in 20170824:

$ comm -23 20170821.ids 20170824.ids | wc -l
     643

20170824 has no sequences that are not in 20170821:

$ comm -13 20170821.ids 20170824.ids | wc -l
       0

and the two have 9111 sequences in common:

$ comm -12 20170821.ids 20170824.ids | wc -l
    9111

These numbers match with the numbers from wc -l above.

So 20170821 is just 20170824 with an extra 643 sequences.

Maybe LC074724 matches, but with a low bit score?

BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000-matches.txt

And....

$ head 20170821-10000-matches.txt
Q   Q   3301
Q   X234A   2616
Q   CS388973.1  2590
Q   DM059402.1  2590
Q   LC074724.1  2581
Q   KJ843188.1  2576
Q   AB697496.1  2576
Q   AB697499.1  2576
Q   AB697503.1  2576
Q   AB697508.1  2576

There it is! And it's in 5th position, just as it should be (based on the bitscores discussed above).

Just to confirm, taking away the -max_target_seqs option gets no match:

$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
$

Three confusingly named BLAST options

I've always found the names and descriptions of the following three BLAST command-line options very confusing.

From the Formatting options section of blastn -help:

  -num_descriptions =0>   Number of database sequences to show one-line descriptions for   Not applicable for outfmt > 4
  Default = `500'    * Incompatible with:  max_target_seqs

  -num_alignments =0>   Number of database sequences to show alignments for   Default = `250'    * Incompatible with:
  max_target_seqs

And from the Restrict search or results section:

  -max_target_seqs =1>   Maximum number of aligned sequences to keep   Not applicable for outfmt <= 4   Default = `500'    *
  Incompatible with:  num_descriptions, num_alignments

From the section names, it seems like only the latter option (-max_target_seqs) would have any effect on the search, and that the former two are just about what is displayed.

But... using either just -num_descriptions or -num_alignments with a big value also gets a match:

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_descriptions 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_alignments 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

So BLAST is apparently finding the LC074724 match even without the high -max_target_seqs option, but it's not displaying it unless there is a high -num_alignments or -num_descriptions value. (Keep in mind, this is the 5th best hit! Not somewhere >500)

Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$

Tl;dr:

In summary, the options

 -evalue
 -max_target_seqs
 -num_alignments
 -num_descriptions

are NOT just filters applied to the results, but instead directly influence the results found.

This means you can End up with a bacterial sequence as BEST hit when using a value x, and with no bacterial sequences at all in your top20 when using value y.

To quote the NCBI user service:

  We don't consider this a bug, but I agree that we should document this
  possibility better. This can happen because limits, including max
  target sequences, are applied in an early ungapped phase of the
  algorithm, as well as later.

And

  The cutoffs are applied during an ungapped extension phase, and some
  alignments may improve in score during gapped extension. The safest
  approach is to relax (increase) the threshold and deal with a larger
  result set.

This was totally unexpected behaviour for me and is as far as I know is nowhere documented in the BLAST help or similar.


What does the "X_OS_IND" column mean? I want to survival analysis using the subset of TCGA LUAD dataset, which identifier are located here.

> head(LUAD_clinicalMatrix[,c("X_EVENT","X_OS_IND")])
  X_EVENT X_OS_IND
1      NA       NA
2       0        0
3       1        1
4       0        0
5       0        0
6       0        0

> all(LUAD_clinicalMatrix$X_EVENT==LUAD_clinicalMatrix$X_OS_UNIT)
[1] FALSE

I want to know which identifier means the sample is alive or dead, and "1" means alive or dead?

According to this link, 1=death 0=censor and null=no data.


How can I only get the species name for fasta sequences from blast results? I am trying to make a phylogenetic tree from sequences obtained with blast. I have files that contains hundreds of fasta protein sequences which are all named like >NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]. But if I use these files to do multiple alignments or phylogenetic trees using either ClustalX or MEGA7, every sequence are named after the first ID in the sequence names (for this example it would be NP_567483.1), and what I want is to have the name of the strain instead (here Pseudomonas aeruginosa PAO1) which is part of the sequence description.

So I don't know whether the softwares I'm using understand that the information between the [ ] is the species and that I just didn't find where to set this as the name or if I have to write some kind of a script to change my fasta files before. I'm really new to this and don't know much about it. Thanks for your help in advance. Julian

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

sed 's/ /_/g'

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

Hope this is helpful! :)


alternatives to PHASTER phage annotation server I would like to use PHASTER (http://phaster.ca) from the University of Alberta, but it's currently down for maintenance, as is its predecessor, PHAST.

What alternatives to PHASTER are there for searching and annotating phage sequences within assembled genomes?

EDIT:

I have seen that PROKKA has some entries mentioning the possibility of using PHASER (currently discarded as an option) but also people using their own HMMs to annotate phage sequences using PROKKA, e.g.:
https://github.com/tseemann/prokka/issues/375

Any ideas?

The answer is to use phastaf, a recently developed software tool developed by Torsten Seemann (the same developer as prokka), which comes with a packaged phaster.ca db.  Get it here


Inspection of gene expression in scRNA-seq data I am running the data preprocessing pipeline for scRNA-seq data presented here.

3.8.6.1 Gene expression
In addition to removing cells with poor quality, it is usually a good idea to exclude genes where we suspect that technical artefacts may have skewed the results. Moreover, inspection of the gene expression profiles may provide insights about how the experimental procedures could be improved.
It is often instructive to consider the number of reads consumed by the top 50 expressed genes.

By analyzing my data I have the following most expressed 50 genes (evaluating the total number of transcripts across cells):
Original

Log-transformed

Since I'll be working on log-transformed data, I assumed that the distributions of the top 50 genes were not excessively skewed. However, as you can see in the first plot, there are some samples (cells) with very high number of transcripts for the top gene. Is it something I should worry about? What other insights can I get from these two figures?
Edit
Original analysis (figure attached) was:

The distributions are relatively flat indicating (but not guaranteeing!) good coverage of the full transcriptome of these cells. However, there are several spike-ins in the top 15 genes which suggests a greater dilution of the spike-ins may be preferrable if the experiment is to be repeated.

In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples prior to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the point the course instructors are trying to get across with that message box.

Eyeballing the results, you should also exclude any ERCC genes as well, as they are synthetic constructs and not from the target samples (although record their counts, as it may be informative for QC/normalisation purposes); they'll have odd distributions because of this.

ERCC sequences are typically added based on the concentration of DNA in a sample, rather than the number of cells. However with single-cell data (i.e. with known cell counts), it'd be possible to do a proper ERCC-based normalisation, assuming an equal amount were added for each cell. I suspect that knowing a bit more about the ERCC addition protocol would help in understanding why their distribution is so different.


MSA (protein) with biopython or something else? I am very new to bioinformatics (and python in general), but I would like to use python to more efficiently analyse enzymes both in terms of structure and functio, using Jupiter notebook. I would like to ask what is the best program/source code for multiple sequence alignments (amino acids) to identify conserved binding sites etc. I see that biopython has a few ways of creating alignments, but I have the impression that it is more focused on nucleotide sequences?

In my personal experience, MUSCLE is the easiest program to use in conjunction with Biopython. Biopython features a command line wrapper for this program, which makes it very easy to use. Make sure to download the appropriate MUSCLE program from drive5 and save it somewhere. E.g., if you are using Jupyter in Linux:
!wget https://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz
!tar -xzvf muscle3.8.31_i86linux64.tar.gz
!cp muscle3.8.31_i86linux64 /usr/local/bin
!chmod 755 /usr/local/bin/muscle3.8.31_i86linux64

Then, you can run MUSCLE like so:
from Bio.Align.Applications import MuscleCommandline
def runMUSCLE(infile, outfile):
    muscle_exe = r"/usr/local/bin/muscle3.8.31_i86linux64" #Here is where we installed MUSCLE
    muscle_cline = MuscleCommandline(muscle_exe,
                                     input=infile,
                                     out=outfile,
                                     clwstrict=True #Output in clustal format (more visually pleasing), otherwise the output is in FASTA. Whichever you need.
                                    )
    muscle_cline()

You don't need to specify that your sequence is amino acidic, however keep in mind that the input file must be in FASTA format.


How to compute LD for pairs of variants with Plink I have a set of SNPs (single nucleotide polymorphisms) { S1, S2, ..., SN } from approximately 200 humans. I wish to determine the linkage disequilibrium (LD) for a defined subset of these SNPs. By this I mean e.g. I would like to compute e.g. LD between S1 and S2, S3, S4, but not the LD between S2 and S3.

Could this calculation be performed in PLINK?

Yes, you can do this with plink.

You can use the flag --ld-snp <variant_id> to calculate LD between your SNP of interest and all other SNPs.
If you have multiple SNPs of interest, and want each of those compared with every SNP in your dataset, use --ld-snps <snp1-snp100,snp102,snp104>
A final option to calculate LD between a specific pair of SNPs is --ld <variant_id1> <variant_id2>

This the relevant plink manual page


List of human transcription factors I'm looking for a list of human transcription factor (TF) genes.

Where can I find a TF database or how can I retrieve this list from other databases?

Lambert et al. (2018, Cell. 172(4):650-665) provide a manually curated list of human TFs at this website.


Difference between copy number neutral reads and active reads In this paper, the authors talk about copy number neutral reads (as reads that equally represent parental chromosomes) and active reads (as reads from only one parent chromosome):

  We reasoned that the reads mapped into a genomic window can be partitioned in two sets: one set includes reads that equally represent parental chromosomes (copy number neutral reads); and the other set contains reads from only one parent chromosome (active reads).

Nonetheless, I have some difficulties understanding these two concepts. Can somebody try to expand on what they wrote?

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

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


Difference between genometools installed with anaconda vs apt-get? I installed genometools with anaconda conda install -c bioconda genometools. Something is installed though apparently no binary called gt. If I type in gt Ubuntu suggests to install genometools via apt. What is the difference between both versions? Why is there no gt binary in the anaconda package?

These appear to be two separate collections of tools. The conda version is a collection of python scripts whereas the apt version is based on a C library. So, as you might expect, the python version does not have a gt binary.

They also have separate GitHub accounts conda genometools vs apt genometools.


How to perform bowtie2 analysis with slurm? I am trying to run my alignment script that works locally, using SBATCH.
Official manual for bowtie2 says I can use samtools for that. However when I run the following script I have an error
line 22: bowtie2: command not found

Below is the script:
#!/bin/bash

#SBATCH --time=05:00:00
#SBATCH --account=def-myaccount
#SBATCH --mem=200000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load samtools/1.9

for FILE in *R1_001.fastq.gz;

do
        BASENAME=${file%%_R1_001.fastq.gz*}

        NUM=`echo $FILE | sed 's/_S.*$//'`

        bowtie2 -x cov_nastya -U $FILE -S ${NUM}exit.sam

done

It seems you do not have bowtie installed or loaded. Use module load bowtie2 together when you load samtools.


Too many requests error when querying NCBI I keep getting a "too many requests error" when querying the NCBI SRA database, even though I'm running less than 10 requests per second, and I have an API key, which supposedly should allow me to run 10 per second.
Here is my (Python) code:
import subprocess
import concurrent.futures
import time

bioprojects = [
    "PRJNA644722",
    "PRJNA644892"
]

all_metadata_generators = []

def fetch_metadata(bioproject):
    return subprocess.run(f"esearch -db sra -query '{bioproject}[bioproject]' | efetch -format runinfo", shell=True, capture_output=True, text=True)

with concurrent.futures.ThreadPoolExecutor() as executor:
    for i in range(0, len(bioprojects), 10):
        all_metadata_generators.append(executor.map(fetch_metadata, bioprojects[i:i+10]))
        time.sleep(1)

for metadata_generator in all_metadata_generators:
    for metadata in metadata_generator:
        print(metadata)

The execution of your EDirect commands may be limited unless you are using an API Key, which you can create from your MyNCBI account. Your API Key can be used in one of several different ways.
In an argument to an EDirect command:
esearch -db nuccore -query 'some query' -api_key 12345

In a UNIX environment:
You can either run this just before you are running some edirect commands or add the following line to your .bashrc file.
export NCBI_API_KEY=12345

In the environment, from a Perl script:
$ENV{NCBI_API_KEY} = "12345"


Grouping a set of genes in a ggplot object I have a boolean matrix of mutated or non mutated genes across some samples

> head(muts)
                  TP53 CDKN2A TSHZ3 ARID1A EPHA3 NOTCH1 ACVR2A ARID1B ARID2
LP6008336-DNA_G01    1      1     1      0     0      0      0      0     0
LP6008460-DNA_D01    1      1     0      1     0      0      0      0     0
LP6005500-DNA_D03    1      1     0      0     0      0      0      0     0
LP6008334-DNA_B02    1      0     1      0     0      0      0      0     0
LP6008334-DNA_A04    1      0     1      0     0      0      0      0     0
LP6005334-DNA_H01    1      0     0      0     1      0      0      1     0

By this code

pw_rt = data.frame(
  gene = c('PIK3CA', 'EGFR', 'ERBB2', 'PTEN'),
  pw = c('RTK'),
  stringsAsFactors = F
)

pw_cc = data.frame(
  gene = c('CDKN2A', 'CCNE1'),
  pw = c('Cell Cycle'),
  stringsAsFactors = F
)

pw = rbind(pw_rt, pw_cc)

pw_col = RColorBrewer::brewer.pal(n = 3, name = 'Set1')[1:2]
names(pw_col) = unique(pw$pw)

require(dplyr)

lmuts = reshape2::melt(muts) %>% as_tibble()
colnames(lmuts) = c('sample', 'gene', 'value')

pw_n = pw$pw
names(pw_n) = pw$gene

Np = nrow(muts)

lmuts %>%
  mutate(gene = paste(gene), PW = pw_n[gene]) %>%
  group_by(gene, PW) %>%
  summarise(N = sum(value == 1)) %>%
  ungroup() %>%
  mutate(
    PW = ifelse(is.na(PW), "None", PW),
    N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
    ) %>%
  ggplot(aes(x = PW, y = gene, fill = PW)) +
  geom_tile() +
  geom_text(aes(label = N)) +
  theme_light() +
  theme(legend.position = 'bottom') +
  guides(fill = guide_legend('Pathway')) +
  scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
  scale_x_discrete(limits = c(names(pw_col), 'None')) +
  labs(
    x = "",
    y = 'Gene',
    title = "Occurrence of mutations"
  )

I am able to get proportion of altered pathway by these genes like

However I want to highlight oncogenes and tumour suppressors within this plot

For instance TP53 is a tumour suppressor and CCNE1 is an oncogene

I don't know how to manipulate the above code to highlight or partition genes to tumour suppressor or oncogene

Can you help?

EDITED

I have add extra column to lmuts like

> head(lmuts)
             sample gene value gene_class
1 LP6005690-DNA_H02 TP53     1        TSG
2 LP2000333-DNA_A01 TP53     1        TSG
3 LP6005409-DNA_D03 TP53     1        TSG
4 LP6008141-DNA_H02 TP53     1        TSG
5 LP6008336-DNA_E02 TP53     1        TSG
6 LP6008269-DNA_B06 TP53     1        TSG
>

but

> lmuts %>%
+     mutate(gene = paste(gene), PW = pw_n[gene]) %>%
+     group_by(gene, PW) %>%
+     summarise(N = sum(value == 1)) %>%
+     ungroup() %>%
+     mutate(
+         PW = ifelse(is.na(PW), "None", PW),
+         N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
+     ) %>%
+     ggplot(aes(x = PW, y = gene, fill = PW)) +
+     geom_tile() +
+     geom_text(aes(label = N)) +
+     theme_light() +
+     theme(legend.position = 'bottom') +
+     guides(fill = guide_legend('Pathway')) +
+     scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
+     scale_x_discrete(limits = c(names(pw_col), 'None')) +
+     labs(
+         x = "",
+         y = 'Gene',
+         title = "Occurrence of mutations"
+     )+facet_wrap(~gene_class, ncol=1)
Error: At least one layer must contain all faceting variables: `gene_class`.
* Plot is missing `gene_class`
* Layer 1 is missing `gene_class`
* Layer 2 is missing `gene_class`
Run `rlang::last_error()` to see where the error occurred.

If you had the classification into oncogene/TSG for each gene in a column (let's call it gene_class) in lmuts, you could split both groups into facets within the plot by adding:

+ facet_wrap(~gene_class, ncol=1)


sva for RNA-Seq data without known phenotype I have been working on RNA-Seq data from two different cohorts, and they show very strong batch effect (~35% variance explained by 1st component in PCA). Since I am trying to do a class discovery from a data set with the subtype of only some samples are known, the only methods I have been using are ComBat and pSVA from SVA package. I would like to compare them with more methods to make the clustering result more trust-worthy.

I found that some papers used sva package, but not citing the paper on pSVA, meaning perhaps there are some other ways to use sva for this type of data that I am not aware of. Is there a method to construct a null model and a full model for sva only by knowing the batch the samples are from?

Any other suggestion of batch effect correction tools for a data set with no known phenotype/subtype are also welcome.

You can just remove the first component from the dataset by setting the first eigenvalue to 0 in the diagonal matrix and then multiplying the SVD matrices.

I am not sure which R code you are using to estimate the components (I have a comparison of PCA code in R on https://github.com/aedin/ODSC_2018/tree/master/PCA_vignette)

In ade4 R package, there is a function reconst() which reconstructs data from a selection of components: https://rdrr.io/rforge/ade4/man/reconst.html

Happy to provide code to do this in other PCA versions if needed
Aedin Culhane


How to colour multiple residues in Pymol? We can use

color yellow, resi 3
color yellow, resi 3-8

to colour one residue or a range of residues. But what if I want to colour residue 3-8 and 10-20, how to write it?

Also, my protein has two chains, light chain (L) with 100 residues and heavy chain (H) with 100 residues. How to colour the residue 10-20, 30-40 of light chain, and the residue 5-10, 20-30 of heavy chain in one go?

Note: as I don't have any working version of pymol at the moment,
I am not able to test the solution properly.
However, an easy way is to first create a selection, and then to color it:

select toBecolored, resi 3-8
color yellow, toBecolored

Regarding the L/H chain part, if they are effectively annotated as L and H, something like the following should work:

select toBeColored, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)
color yellow, toBeColored

Combining both commands is probably possible as:

color yellow, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)

Pymol allows synonyms (eg. resi and i.  or chain and c.).

The following pages may be also helpful:

Pymol wiki: help on Selection
SourceForge

And the two following links, if you want to get a bit deeper in proper programming with python for pymol (as they can provide you with inspiration):

Pymol wiki: plugins
Pymol scripts


What will be an appropriate mathematical distribution for SNP data? I found that several papers describe SNPs as a binomial distribution with the probability of "success" equals to minor allele frequency.

However, in my experiments, when I generate SNP array following this distribution, the simulation results behave very differently from the SNP array generated following the random mating procedures. I feel like there is some missing information about the description of a binomial distribution.

I wonder what will be a better mathematical description of the SNP array data? Besides a direct answer, any thoughts, or any suggestion of papers are also highly welcome.

EDIT:

Thanks for these suggestions, after reading relevant materials, I think what I observe is about the Hardy-Weinberg equilibrium affected by evolutionary influences like mutation and migration. I will continue to look and list findings here if I find anything interesting. Please feel free to suggest some mathematical descriptions.

Also thanks for the suggestion of tools. They look very helpful, but I really hope to understand some mathematics behind it.

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

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


CIBERSORT runtime error I was running CIBERSORT web and saw this error:

Exception in thread "Thread-0" java.lang.ArrayIndexOutOfBoundsException: 0
at launch.QN.doQN(QN.java:56)
at launch.Cibersort.Run_CIBERSORT(Cibersort.java:611)
at launch.Cibersort.run_all(Cibersort.java:179)
at launch.Cibersort$runCIBERSORT.run(Cibersort.java:1585)
at java.lang.Thread.run(Thread.java:745)
311

Does anyone know why?

OK I found the bug. CIBERSORT seems to require tsv files instead of csv files. Plus it doesn't seem to like empty columns


What is the difference between the entity_src_genCategory tag and the entity_src_natCategory tag in pdb-xmls? I want to read out the taxonomy id of the source organism of a pdb entry using the xml-file of the entry that is provided under http://files.rcsb.org/view/{pdb_accession}-noatom.xml.

I have so far found the taxonomy id in the tag PDBx:entity_src_genCategory (example entry 1UKV) but also in the tag PDBx:entity_src_natCategory (example entry 2BWI).

XML for 1UKV with PDBx:entity_src_genCategory

<PDBx:entity_src_genCategory>
  <PDBx:entity_src_gen entity_id="1" pdbx_src_id="1">
     <PDBx:gene_src_common_name>baker&apos;s yeast</PDBx:gene_src_common_name>
     <PDBx:gene_src_genus>Saccharomyces</PDBx:gene_src_genus>
     <PDBx:host_org_genus>Escherichia</PDBx:host_org_genus>
     <PDBx:host_org_species>Escherichia coli</PDBx:host_org_species>
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_gene_src_ncbi_taxonomy_id>4932</PDBx:pdbx_gene_src_ncbi_taxonomy_id>
     <PDBx:pdbx_gene_src_scientific_name>Saccharomyces cerevisiae</PDBx:pdbx_gene_src_scientific_name>
     <PDBx:pdbx_host_org_ncbi_taxonomy_id>469008</PDBx:pdbx_host_org_ncbi_taxonomy_id>
     <PDBx:pdbx_host_org_scientific_name>Escherichia coli BL21(DE3)</PDBx:pdbx_host_org_scientific_name>
     <PDBx:pdbx_host_org_strain>BL21(DE3)</PDBx:pdbx_host_org_strain>
     <PDBx:pdbx_host_org_vector_type>PLASMID</PDBx:pdbx_host_org_vector_type>
     <PDBx:plasmid_name>pET19tev</PDBx:plasmid_name>
  </PDBx:entity_src_gen>

XML for 2BWI with PDBx:entity_src_natCategory

<PDBx:entity_src_natCategory>
  <PDBx:entity_src_nat entity_id="1" pdbx_src_id="1">
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_ncbi_taxonomy_id>223</PDBx:pdbx_ncbi_taxonomy_id>
     <PDBx:pdbx_organism_scientific>ACHROMOBACTER CYCLOCLASTES</PDBx:pdbx_organism_scientific>
  </PDBx:entity_src_nat>

What is the difference between them? Are there other tags that contain information about the source organism of the entry? Can they occur together in one file? If so, is one of them more valid?

src_nat = isolated directly from a natural tissue

src_gen = the source was genetically manipulated

The files that you use (PDBML/XML) are derived from the primary PDB data format which is PDBx/mmCIF.
The names in XML seem to correspond to mmCIF names, so you can read documentation on http://mmcif.wwpdb.org/, and here are direct links:

category entity_src_nat
category entity_src_gen

Can they occur together in one file? I suppose they can for complexes: one "entity" could be natural and another genetically manipulated.


Bootstrap analysis - why it is called bootstrap? Bootstrap analysis, bootstrapping etc are quite common jargons of bioinformatics and phylogenetics.

However, it is not very much clear to us, what exactly being meant by "a boot's straps".

Does it means a "comparison"? (Such as we hold the two straps of a boot (shoe) in "two hands"; may be an analogy to comparing 2 things (say left-'hand' side and right-'hand' side of an algebraic proof-task)) or it is something about sequence alignment (since we 'tie' a knot in boot (shoe) so that the 2 ends of the strap remain attached to each-other).

Some of my classmates even assumed it may be a process during computer booting (but I could not agree with that).

I've searched internet, but could not found anything very helpful. This wikipedia page (permalink) about bootstrapping, tells,

   Tall boots may have a tab, loop or handle at the top known as a bootstrap, allowing one to use fingers or a boot hook tool to help pulling the boots on. The saying "to pull oneself up by one's bootstraps" was already in use during the 19th century as an example of an impossible task. The idiom dates at least to 1834, when it appeared in the Workingman's Advocate: "It is conjectured that Mr. Murphee will now be enabled to hand himself over the Cumberland river or a barn yard fence by the straps of his boots." In 1860 it appeared in a comment on philosophy of mind: "The attempt of the mind to analyze itself an effort analogous to one who would lift himself by his own bootstraps".  Bootstrap as a metaphor, meaning to better oneself by one's own unaided efforts, was in use in 1922. This metaphor spawned additional metaphors for a series of self-sustaining processes that proceed without external help.

  The term is sometimes attributed to a story in Rudolf Erich Raspe's The Surprising Adventures of Baron Munchausen, but in that story Baron Munchausen pulls himself (and his horse) out of a swamp by his hair (specifically, his pigtail), not by his bootstraps – and no explicit reference to bootstraps has been found elsewhere in the various versions of the Munchausen tales.

(In the quotation, highlighted portion seemed to tried to tell a meaning), but it wasn't also very much helpful. Could not correlate with bioinformatics. (Btw the same wikipedia article mentions computer booting).

Wikipedia article about Bootstrapping (statistics) (permalink) tells:

  The bootstrap was published by Bradley Efron in "Bootstrap methods: another look at the jackknife" (1979), inspired by earlier work on the jackknife.Improved estimates of the variance were developed later. A Bayesian extension was developed in 1981. The bias-corrected and accelerated (BCa) bootstrap was developed by Efron in 1987, and the ABC procedure in 1992.

Still that is not very helpful.

However, the textbook, Introduction to Bioinformatics/ Arthur M. Lesk/ 3rd Edition/ Oxford / (Low Price edition) ; in its chapter 5 (Alignments and phylogenetic trees), section The problem of varying rates of evolution -> computational consideration (page 296); a clearer definition has been given:

  Formal statistical tests, involving rerunning the calculation on subsets of the original data, are known as jackknifing and bootstrapping.

  jackknifing is calculation with data sets sampled randomly from the original data. For phylogeny calculations from from multiple sequence alignments, select different subsets of the position in the alignment, and rerun the calculation. Finding that each subset gives the same phylogenetic tree lends it credibility. If each subset gives a different tree, none of them are trustworthy.
  Bootstrapping is  similar to jackknifing except that the position chosen at random may include multiple copies of the same position, to form data sets of the same size as original, to preserve statistical properties of data sampling.

More clear definition of "Boot-strapping"; but it does not explain, where the relation exists with the boot (shoe) and its straps.

The method itself has nothing really to do with boots or straps, just as the jack-knife method also has nothing to do with knives. In the case of bootstrapping, the goal is to determine the accuracy of an estimate from random subsets of a population. Normally estimating something like the variance of the mean requires multiple independent samples, but bootstrapping allows you to perform estimates from a single population. This is essentially allowing the estimate to "pull itself up by its bootstraps", which is the only reasonable source for the term. See also the more general wikipedia article on bootstrapping outside of statistics.


freebayes segmentation fault I have a draft, Illumina-based, genome and bunch of reseq samples. I generated bam files using bowtie2 and marked duplicates using picard.
Now, when I tried to call SNPs using freebayes, I got cryptic Segmentation fault (core dumped) even when I tried to run it in verbose debugging mode (-dd). Do you have an idea where might be a problem?
freebayes -f data/Afus1/genome.fa -b reseq1.sorted.rmdup.bam reseq2.sorted.rmdup.bam -dd
  loading fasta reference data/Afus1/genome.fa
                                                                         Opening 2 BAM format alignment input files

     done
Number of ref seqs: 821466
Number of target regions: 0
no sample list file given, reading sample names from bam file
found 2 samples in BAM file
Segmentation fault (core dumped)
```

Oh you silly sausage!
This obviously means that you are running out of stack memory. Apparently...

freebayes is creating a single string containing the entire vcf header, which gets quite big when the reference has a lot of contigs with long names, causing stack space depletion.

Try increasing the maximum stack size (in this example to 128M) and running your command again.
ulimit -S -s 131072


Determining significance of a variable in a glm model How would one determine the significance of a variable in a glm model?
If I, for example, have a dataframe like seen below, how would I determine if the origin of the sample has a significant effect on the value? (this is the number of enzymes capable of degrading the substrate f that matters)
Substrate    variable value origin
cellulose       M09    8    free
mannan          M12    2    free
glycogen        M65    2    free
chitin          M87    4    free
cellulose       M90    2    isolate
manan           M78    1    isolate
glycogen        M21    4    isolate
chitin          M21    1    isolate

So far I have tried:
mcomp = glm.nb(value ~ origin, data = my_data)

summary(mcomp)
Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.9625  -0.9047  -0.9047   0.1212   3.5232

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.01657    0.06571  -0.252  0.80097
originisolate -0.21911    0.08180  -2.679  0.00739 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.3418) family taken to be 1)

    Null deviance: 2053.5  on 2679  degrees of freedom
Residual deviance: 2046.3  on 2678  degrees of freedom
AIC: 6517.5

Number of Fisher Scoring iterations: 1

              Theta:  0.3418
          Std. Err.:  0.0186

 2 x log-likelihood:  -6511.4590

So free becomes the intercept and then isolate if significantly different from that. Does this mean Origin has a significant effect on the value?
Would the better approach be to do the following?:
mcomp = glm.nb(value ~ origin + Substrate, data = comb_data)
summary(aov(mcomp))
              Df Sum Sq Mean Sq F value Pr(>F)
origin         1     23   22.55   6.612 0.0102 *
Substrate     44   1445   32.84   9.631 <2e-16 ***
Residuals   2634   8981    3.41
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

This shows me that origin and substrate have an effect on value if I understand correctly?

There is no better method, it's a matter of what you want to test or what is your question.
Using the anova() or aov(), test the terms collectively. For example, in your example with Substrate, the null hypothesis is that the coefficients are all zero, meaning cellulose =0, mannan =0 , ....
If the question is, "do the isolate samples have a higher value than origin samples?", then you can use your first model, where free is set as the reference and you test whether the effect of isolate is non-zero. Likewise you can do this for substrate and set of them as your reference. You can also do other pairwise comparisons using this model.
If the question is, "does origin have a significant effect on value, after controlling for substrate?", then you can use your second model.


Pre-filtering genes for Principal Component Analysis I have a raw counts data-set of 20,502 genes and 137 samples. I want to find out Principal Components which best explain variation between samples in different stages of tumor.

I am new to Machine learning and would like some help in selecting a criteria to reduce the number of genes before doing a PCA.

I have tried doing a DESeq2 analysis on the data-set but I am confused as to what to do next.

You can filter out genes that do not show variation across samples, they would not differentiate samples anyway. For the specifics on how to do so, please see the Subpopulation Analysis section of a nice single cell RNA-seq workflow from the Kharchenko Lab. I believe there won't be fundamental differences when applying this approach to bulk RNA-seq.


How to represent pathways as mathematical graphs? I would like to construct a mathematical graph: $G = (E,V, w(E))$ a graph with vertices, edges with weight on them based on a KEGG pathway.

Is there a way to do that with existing tools? If not I will start to implement a Python module for that. My main question is what is the mathematical graph representation of pathways? Assuming the possible node types and the different connections?

I would like to transform a pathway to a graph because then one can run min-cut-max-flow algorithms and stuff like that as part of my pathway analysis.

Pathways are generally directed graphs, which can be implemented in python using the networkx package. Note that there will be cycles in graphs.

Please note that there are large numbers of pathways analysis packages that already exist (both python and R packages but also just webservers), so I strongly recommend searching to see if something already exists before wasting time and implementing it yourself.


Do people today with a disease caused by single gene mutation that is inheritable - all share a common ancestor? First question here so sorry if it is about both genetics and evolution.

So, if I get a disease which can be clearly identified on the gene level and all my descendants inherit it, after 500 years do all people with the disease share the common ancestor (given that it's only inheritable)?

I mean, my clear guess first would be "yes" but I really want to check it since I haven't read it anywhere.

Feel free to post materials on this too.
Thanks!

That can happen, but it's more likely that the disease will be present in more than one family. Diseases with a single, small genetic cause will be more likely to flare up elsewhere (assuming the disease isn't too bad, and it wouldn't be if it's been in the family for 500 years).

However, what you will find is that in your family that disease is almost always present together with some other inherited genetic feature that is nearby on the genome (i.e. an inherited haplotype block). This association is less likely to be shared with other families.


How can I assign a genomic region into a window using R? I have asked this question at BioStars but have not gotten any suggestions so far. So, I am asking the community here, perhaps someone here has an idea how I can solve my problem.

I am dealing with a bunch of inversions and running some analysis on each of them. One of them is located on chromosome one, the starting position is 6199548 and ending position is 9699260. So the length of this variant is 3499712 bp. I want to generate non-overlapping windows based on the length of these focal positions (3499712) and assign the focal position to one of these windows that overlap with the starting and ending position of my focal position. I am achieving this by using the GenomicRanges::tileGenome function in R:

seqlength <-c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)
bin_length<-c(3499712)
 bins<-tileGenome(seqlength,tilewidth=bin_length,cut.last.tile.in.chrom
        = T)
            #bins_u<-unlist(bins)
            write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names
        = FALSE, row.names = TRUE) }

This generates non-overlapping 3499712 bps windows across the genome.

window_id chromosome start end length
"1" "chr1" 1 3499712 3499712
"2" "chr1" 3499713 6999424 3499712
"3" "chr1" 6999425 10499136 3499712
"4" "chr1" 10499137 13998848 3499712
"5" "chr1" 13998849 17498560 3499712

However, I have a major problem here that I have not been able to sort out. I want one of these windows to overlaps exactly with the starting and ending position of my focal position (6199548:9699260). Only the focal position not any other position outside this region. This means 2 windows (one before and one after the focal position) will have a different length which my analysis is fine with. I just do not know how to do it?

This is my desired output:

window_id chromosome start end length
"1" "chr1" 1 3499712 3499712 "*"
"2" "chr1" 3499713 6199547 2699834 "*"
"3" "chr1" 6199548  9699260 3499712 "*"
"4" "chr1" 9699261 13998848 4299587 "*"
"5" "chr1" 13998849 17498560 3499712 "*"

As you see, the the second window is now smaller, the third window accommodates the focal position and the forth window is a little larger than the bin size. Does anyone have any idea how can I do this task in R?

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

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

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

bin_length <- width(site)

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

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

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

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

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

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

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

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

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


calculating nucleotide frequency per column I have some sequences shown below
CAGGTAGCC
CCGGTCAGA
AGGGTTTGA
TTGGTGAGG
CAAGTATGA
ACTGTATGC
CTGGTAACC
TATGTACTG
GCTGTGAGA
CAGGTGGGC
TCAGTGAGA
GGGGTGAGT
TGGGTATGT
GAGGTGAGA
CAGGTGGAG

Each line has 9 nucleotides.
Consider it to be 9 columns.I want to calculate the nucleotide frequency of each nucleotide for each of the 9 columns.
For example 1st column will have these bases C,C,A,T,C,A etc
Out put should be something like this
A   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
C   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
G   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
T   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36

Note, I just made up the numbers to show you the output file format


population admixture I am new to Admixture and I would like to use fastNGSadmix to look for admixture in one individual (modern). Admixture is supposed to look for proportions of ancestral populations in a sample (either a population or an individual sample). for that I will need to prepare my reference population set to compare with my actual sample which is an individual form region X for instance.

Now I am thinking whether I have to add the X population to my reference set or not? If yes, wouldn't it introduce bias? Like my individual will always show the X population as the first admixing population?

Disclaimer, I am not an expert on admixture, answer is provided without any warranty.

It depends what you want to figure out. Why do you want to categorize the ancestry of the individual? If you analyze this origin of population X you should make a reference without population X and analyze all the individuals from the population with all the surrounding populations in the reference. However if you are interested in origin of a single individual, than the population X got to be in reference as well.

If your result will be that 90% of SNPs of your individual correspond to your population X, then it means that's a normal X fellow. Why it should be a problem?

-- edit --

I just found this preprint that highlights potential pitfalls in interpretation of STRUCTURE/ADMIXTURE results.


how to merge more than two sample in Seurat? I would like to merge more than two sample in the Seurat, and the mergeseurat can only merge two sample. So what should I do now. The screenshot is my script.

You have fed arguments to the MergeSeurat() function that it does not expect. In terms of objects, MergeSeurat() accepts only 2 arguments object1 and object2, please do ?MergeSeurat and see yourself.

Seurat is probably the best documented single cell package, your exact question has its own dedicated vignette: https://satijalab.org/seurat/v3.0/merge_vignette.html

As the vignette above is for Seurat v3 and in case you would like to use Seurat v2:
You can create a named list of Seurat objects to be merged and then use reduce() of the purrr package.

# create a named list of seurat objects to be merged
# code not shown

# optional but probably a good idea
# rename cells using object names as prefix
for (i in names(seurat_object_list)) {
  seurat_object_list[[i]] <- RenameCells(seurat_object_list[[i]],
                                         add.cell.id = i)
}

# merge all the objects in the list
merged_combined <- reduce(seurat_object_list,
                          MergeSeurat,
                          do.normalize = FALSE)


RNASeq analysis using featureCount and EdgeR I am using a pipeline (bam -> featurecount-> EdgeR) to do some RNASeq analysis of several groups and sub-groups. For example, I have the following dataset with two types (T1 and T2) and T1 has three groups G1, G2 and G3, T2 has G5 and G6 groups. I would like to find TPM numbers, fold-change and FDR for all the genes and Groups.
Type Group SRA_Run#
T1   G1.   <some SRA ids e.g. SRR123456, SRR112233>
T1   G2.   <some SRA>
T1.  G3.   <some SRA>
T2.  G5.   <some SRA>
T2.  G6.   <some SRA>

I have individual BAM files generated for each SRA Run. I also individual read counts from featureCount for each SRA Run. I was not sure if I need to merge individual read counts for each group or run featureCount using all bam files of a group at once. Does it make any difference in  the analysis that follows.
I also would like to know if I can run EdgeR at once for all the types and groups at once for doing differential gene expression for each group of Type 1 (T1) against every other group of Type 2 (T2). I am bit stuck at this point how to proceed. Any help would be appreciated.

It makes no difference if you process the BAM files one at a time with featureCounts or all together, except that it changes how you have to read the files into R.
You can supply edgeR with lists of contrasts to have it compute fold-changes and p-values for. Please have a look at the edgeR user guide for examples.


How do I add external labels and connecting lines to a Ternary Plot in R? I have plotted a triangular plot,that three dimensions of which represent three herbal strategies. One dimension represents the strategy of C (competitive plant), the second dimension “S” (stress tolerant plants) and the third dimension”R” (ruderal plants), as you can see in the image below,
the points on it represent the plant species. Now I am trying to write name of  the plant species outside the triangle and connect with the arrow to the corresponding point. How do I do this work?
Many thanks in advance as it would help me a lot in my work.

TernaryPlot()
#Plot two stylised plots side by side, and plot data
par(mfrow=c(1, 1), mar=rep(0.3, 4))
TernaryPlot(atip='C%', btip='R%', ctip='S%',
            point='UP', lab.cex=0.8, grid.minor.lines=0,
            grid.lty='solid', col='#FFFFFF', grid.col='GREY',
            axis.col=rgb(0.1, 0.1, 0.1), ticks.col=rgb(0.1, 0.1, 0.1),
            padding=0.08)
data_points <- list("Bromus dantonia" = c(47,   59, 149),
                    "Calamagrosis psoudo phragmatis" = c(90,    102,    63),
                    "Carex diluta" = c(109, 64, 82),
                    "Carex divisa" = c(96,  99, 59),
                    "Carex pseudocyperus" = c(130,  71, 54),
                    "Carex stenophylla" = c(97, 98, 59),
                    "Catabrosa aquatica" = c(100,   5,  150),
                    "Centaurea iberica" = c(124,    85, 46),
                    "Cirsium hygrophilum" = c(158,  42, 55),
                    "Cladium mariscus" = c(159, 96, 0),
                    "cod2" = c(54,  82, 119),
                    "Cynodon dactylon" = c(121, 54, 80),
                    "Eleocharis palustri" = c(124,  100,    31),
                    "Epilobium parviflorum" = c(67, 80, 107),
                    "Eromopoa persica" = c(83,  15, 157),
                    "Funaria cf.microstoma" = c(8,  0,  247),
                    "Glaux maritime" = c(4, 196,    55),
                    "Hordeum brevisubulatum" = c(76,    70, 109),
                    "Hordeum glaucum" = c(40,   79, 136),
                    "Inula britannica" = c(95,  108,    51),
                    "Juncus articulatus" = c(107,   79, 69),
                    "Blysmus compressus" = c(81,    127,    47),
                    "Juncusinflexus"= c(149,    106,    0),
                    "Medicago polymorpha" = c(60,   86, 109),
                    "Mentha spicata" = c(150,   23, 82),
                    "Ononis spinosa" = c(66,    112,    77),
                    "Phragmites australis" = c(234, 0,  21),
                    "Plantago amplexicaulis" = c(108,   83, 64),
                    "Poa trivialis" = c(90, 28, 138),
                    "Polygonum paronychioides" = c(20,  12, 223),
                    "Potentila reptans" = c(106,    41, 108),
                    "Potentilla anserina" = c(105,  58, 91),
                    "Ranunculus grandiflorus" = c(129,  25, 101),
                    "Schoenus nigricans" = c(143,   91, 21),
                    "Setaria viridis" = c(10,   7,  238),
                    "Sonchus oleraceus" = c(178,    0,  77),
                    "Taraxacum officinale" = c(117, 28, 110),
                    "Trifolium repens" = c(94,  4,  157),
                    "Triglochin martima" = c(63,    96, 95),
                    "Veronica anagallis-aquatica" = c(55,   37, 163)
)
AddToTernary(points, data_points, pch=21, cex=1.2,
             bg=vapply(data_points,
                       function (x) rgb(x[1], x[2], x[3], 128, maxColorValue=255),
                       character(1))
)

AddToTernary(text, data_points, names(data_points), cex=0.8, font=1)

Hmmm I'm not familiar with AddtoTernary function, I can see that you've placed a function within the "AddtoTernary" function. The traditional ggplot2 way of doing this would simply be controlled by the theme

ggtern(data=x,aes(x1,x2,x3)) +
  geom_point(fill="red",shape=21,size=3) +
  theme_tern_bw()

This will give an arrow for each axis, but additional code is required for different colours (I'd have to look up my past work).

I think what you might want is a colour key? In this case ggtern will have a "key" function because ggplot2 has one.Here, your colouring is based on species not on a colour gradient as to how far to the corner of a graph you are.

Alternatively you can easily put names besides each point in ggplot2, but again this is traditional style code.

Also you can add different shapes to the graph, which might be useful instead of colours.

Arrows and names together peronsally I've not seen on a ggplot2 graphic and I'd have not idea how to do it.

FYI I don't think we need all the data just a few points fine, although it shows why you prefer arrows - personally I would use colour coding. In "traditional" ggplot2 coding the graphing is built up in layers, so if additional functionality is required it is simply appended to the exisiting command. Thus it is quite easy to make modifications.


Inflated p-values in quantitative trait analysis I am performing a quantitative trait association between the expression of one gene and ~400,000 methylation values. First, both variables are rank inverse normal transformed, adjusted for confounders (technical and for family structure) in a mixed-effects model and we use the residuals for the association. Then we do a linear regression between the residuals of expression and each methylation site, adding other four confounders in the model (age, BMI, etc.). The sample size is of ~500. For the association we are using MatrixEQTL.

When we check the p-values, they are clearly inflated. Is it something that I have to be concerned about? Did I probably miss something in the analysis? Thanks in advance.

Permutation as suggested by @StupidWolf's comment is essential to understand what's going on. If permutation makes this pattern go away, then you have a problem with your model specification, there's something uncorrected.

If your data are weird, well, that's just how they are. But this argues to me that something else is going on confounding your associations. Cryptic population structure in the data would lead to this kind of pattern, e.g. here. That is the first thing that most reviewers would say, I'd guess.

Some random thoughts:

What does the distribution of your methylation data look like after normalization?
Your observed distribution is basically never on the expected line, which argues strongly for confounding.
How are you adjusting for non-independence of methylation types? Family-wise correction is likely not specific enough. You may have to estimate an identity-by-state SNP similarity matrix for your population and fit that as part of your model, or at least fit 10-20 principal components to account for population structure.
How are you representing your methylation status? If it's binary you might look at this paper.


Order of batch effects removal, data imputation and library size normalization in scRNA-seq data I am preprocessing scRNA-seq data. What is the best practice in use to run both ComBat for batch effects removal, data imputation (to mitigate dropout) and library size normalization?

I thought that library size should be run first, since it is per-cell normalization, then ComBat batch effects removal. On the original paper - Johnson et al. (2007) - it is stated that:

  We assume that the data have been normalized and expression values have been estimated for all genes and samples.

However, I want to apply it to scRNA-seq data. Does this statement still hold? Additionally, I plan to apply imputation (e.g. with MAGIC) in the end. Is there any problem you can spot?

Update

I attach the PCA regarding an example Mus Musculus dataset in which different colors represent different mice. It seems clear to me that the first two principal components are affected by batches (mouse id).

Update 2

I rerun the PCA on raw counts data (the first PCA was on log-transformed data) and I obtain a different description of the dataset, in which batch effects seem not to be prevalent.

MAGIC assumes input data has been both library-size normalized, and either log or sqrt transformed prior to imputation (see also: MAGIC tutorial). Additionally, any graph-based methods (MAGIC, PHATE, t-SNE, UMAP, spectral clustering, Louvain, etc etc) will give flawed results if your data contains a batch effect, since the neighbourhood graph would reflect that structure of your batch effect, and worse, imputation would further reinforce this batch effect.

Thus I would recommend the following pipeline:

Library-size normalization
Square root (or log) transform
Batch effect removal
Imputation

Regarding your update, the reason you don't see the batch effect in the raw counts data is simply that the batch effect is not visible in the most highly expressed genes. Prior to transformation, the principal source of variation in your data is simply the expression of the most highly expressed genes - this is essentially masking the source of the batch effect, not removing it. I recommend never working with raw molecule counts in scRNAseq, as the raw counts data hides much of the heterogenity in your dataset, which is precisely what you are looking in doing single-cell RNA-seq.