Create a Grouped Barchart I’m trying to create a Group Chart like this:

I have 7 classes for each of 4 species, not 5 like in the diagram. I want to build i diagram to show the distributions of the classes in each species

Here’s example provided online:

  library(ggplot2)

    # create a dataset
    specie <- c(rep("sorgho" , 3) , rep("poacee" , 3) , rep("banana" , 3) , rep("triticum" , 3) )
    condition <- rep(c("normal" , "stress" , "Nitrogen") , 4)
    value <- abs(rnorm(12 , 0 , 15))
    data <- data.frame(specie,condition,value)
    # Grouped
    ggplot(data, aes(fill=condition, y=value, x=specie)) +
        geom_bar(position="dodge", stat="identity")

------------------------------------------
    #Here’s code:

    > Species <- c(rep("Cames1" , 5) , rep("Came2" , 5) , rep("Cams3" , 5) , rep("cams4" , 5) )
    > Class <- rep(c("EVV" , "MaLR" , "psy" , "Nargo" , "Pao" , "Unknown" , "Undefiend") , 8)
    > value <- abs(rnorm(20 , 0 , 15 ))
    data <- read.table(text="Class Cames1 Came2 Cams3 cams4
    + EVV 2194 5435 2706 2902
    + MaLR 503 108 184 231
    + psy 284 187 413 295
    + Nargo 274 1977 413 2195
    + Pao 66 29 78 97
    + Unknown 286 283 273 377
    + Undefined 76 77 80 95")
   > ggplot(data, aes(fill=Class, y=7, x=Species)) + geom_bar(position="dodge", stat="identity")

Error: Aesthetics must be either length 1 or the same as the data (7): x, fill

i am able to build basics chart, but to group them all into one like the example provided challenging for me , could someone help out ?

I have noticed a few errors:

1) You define a bunch of variables (Species, Class, ...) but then instead of creating your data frame data with these you create it by reading from text. I don't think this is good practice, you should do so like in your example code.

2) With the way you define, your data frame does not contain columns like Class or Species but columns V1 to V5, you should use argument header = TRUE in your read.table() call.

3) Even when you read your data as in step 2, you would need to transform your data table into "long format", please see below.

4) In your ggplot() call, you assign 7 to your y aesthetic, this does not make sense. You should either provide a vector of length equal to your other aesthetics (x and fill) or better a column name of values of interest.

# This is your code minus the redundant variables that you have created but are
# not used

library(tidyr)
library(ggplot2)

# Pay attention to the header argument
data <- read.table(
  text = "Class Cames1 Came2 Cams3 cams4
     EVV 2194 5435 2706 2902
     MaLR 503 108 184 231
     psy 284 187 413 295
     Nargo 274 1977 413 2195
     Pao 66 29 78 97
     Unknown 286 283 273 377
     Undefined 76 77 80 95",
  header = TRUE
)

# data in wide format
> data
      Class Cames1 Came2 Cams3 cams4
1       EVV   2194  5435  2706  2902
2      MaLR    503   108   184   231
3       psy    284   187   413   295
4     Nargo    274  1977   413  2195
5       Pao     66    29    78    97
6   Unknown    286   283   273   377
7 Undefined     76    77    80    95

# Making the data wider: Keeping Class column as it is but condensing the data
# within the rest of the columns into two newly created columns, Species_ and
# Value_

data_long <- gather(data,
                    key = "Species_",
                    value = "Value_",
                    -Class)

# First few rows of your data in long format
> head(data_long)
    Class Species_ Value_
1     EVV   Cames1   2194
2    MaLR   Cames1    503
3     psy   Cames1    284
4   Nargo   Cames1    274
5     Pao   Cames1     66
6 Unknown   Cames1    286

# x and y aesthetics are the newly created columns Species_ and Value_
ggplot(data_long, aes(fill=Class, y=Value_, x=Species_)) + geom_bar(position="dodge", stat="identity")


How can I calculate coverage at single bases using a bam file? I'm looking for a way to input a vcf or bed file (with specific base positions) and a bam file, and get the coverage at each base position (ie single base bins) using the bam file. I also want the strand information so ideal output would be something like:

chr base_position total_coverage fwd_coverage rev_coverage

I've tried bamCoverage and bedtools Coverage but with bamCoverage, I can't figure out how to get single-base pair resolution and with bedtools Coverage, I can't figure out how to get the strand coverage information.

Note: not yet tested, so there may be some additional fiddling with command line options needed

The per-base depth can be obtained from samtools depth (-a includes zero-coverage positions):

samtools depth -a in1.bam > depth_in1_both.tsv

To split this by forward and reverse, you can use an initial pipe through samtools view to exclude or include reverse-complement mappings:

samtools view -b -F 0x10 in1.bam | samtools depth -a > depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a > depth_in1_rev.tsv

To filter for particular regions, this can be pre-piped through bedtools intersect:

bedtools intersect -a in1.bam -b regions.bed | \
  samtools view -b -F 0x10 in1.bam | \
  samtools depth -a > depth_in1_regions_fwd.tsv

Combining both, forward, and reverse is a bit trickier. I'd recommend doing it in a more manual fashion in R. But if you want to keep to the command line...

If you are absolutely sure that the 'depth' output will appear in the same order in all files (it might be worth using samtools depth -aa), then the lines can be pasted, with anonymous pipes, with unnecessary fields filtered out using cut:

paste <(samtools depth -aa in1.bam) \
      <(samtools view -b -F 0x10 in1.bam | samtools depth -aa) \
      <(samtools view -b -f 0x10 in1.bam | samtools depth -aa) | \
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

Or equivalently, if using intermediate files that have been generated already:

paste depth_in1_both.tsv depth_in1_fwd.tsv depth_in1_rev.tsv |
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

It would also be possible to use join to combine the files (in cases where file output is not consistent), but the first two fields of each file would then need to be combined, because 'join' only joins on the first field... but if you need to get that complex, it'd be better to use something other than pipes on the command line.


Griding ggplot object in r I have this data frame

  head(d)

# A tibble: 6 x 3
  gene  variable       value
1 CCT5  Res_immune     10.6
2 CCT5  Non_res_immune 10.5
3 CDK1  Res_immune      9.25
4 CXCL5 Non_res_immune  9.82
5 GBP1  Res_immune      9.04
6 GBP1  Non_res_immune  9.43
>

And this code

library(ggplot2)
library(ggrepel)
ggplot(d, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)

The genes in first column are involved in cytokine signaling pathway in two groups of patients

This is my plot

I have another data frame for genes in IL17 signaling in the sample patients like

head(d)

# A tibble: 6 x 3
  gene  variable       value
  <chr> <chr>          <dbl>
1 CDK6  Res_immune     17.6
2 IL6  Non_res_immune 14.5
3 CXCL8  Res_immune      1.25

How I can put plots from both data in a grade like below where the name of pathways shown instead of DEL, INS, SNP and total like below picture?

What you are looking for is facetting. Using this key word on any search engines you will find dozen of answers describing its use: https://ggplot2.tidyverse.org/reference/facet_grid.html and https://ggplot2.tidyverse.org/reference/facet_wrap.html

To prepare the faceting in ggplot2, a possible way is to bind both dataframes together by first specifying which group they belong such as:

cytokine$Group = "Cytokine"
il17$Group = "IL17"
DF <- rbind(cytokine,il17)
DF

    gene       variable value    Group
1:  CCT5     Res_immune 10.60 Cytokine
2:  CCT5 Non_res_immune 10.50 Cytokine
3:  CDK1     Res_immune  9.25 Cytokine
4: CXCL5 Non_res_immune  9.82 Cytokine
5:  GBP1     Res_immune  9.04 Cytokine
6:  GBP1 Non_res_immune  9.43 Cytokine
7:  CDK6     Res_immune 17.60     IL17
8:   IL6 Non_res_immune 14.50     IL17
9: CXCL8     Res_immune  1.25     IL17

Then, use facet_grid function of ggplot2:

library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_grid(Group~., scales = "free")

Alternatively, you can use facet_grid (or facet_wrap) to have your plot in a single row:

library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_wrap(.~Group, scales = "free")

Reproducible data example

structure(list(gene = c("CCT5", "CCT5", "CDK1", "CXCL5", "GBP1",
"GBP1"), variable = c("Res_immune", "Non_res_immune", "Res_immune",
"Non_res_immune", "Res_immune", "Non_res_immune"), value = c(10.6,
10.5, 9.25, 9.82, 9.04, 9.43), Group = c("Cytokine", "Cytokine",
"Cytokine", "Cytokine", "Cytokine", "Cytokine")), row.names = c(NA,
-6L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)

and

structure(list(gene = c("CDK6", "IL6", "CXCL8"), variable = c("Res_immune",
"Non_res_immune", "Res_immune"), value = c(17.6, 14.5, 1.25),
    Group = c("IL17", "IL17", "IL17")), row.names = c(NA, -3L
), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)


coexnet package for network co-expression analysis, is looking for a filelist.txt that doesn't exist On my way to create a co-expression matrix for a microarray experiment, using the coexnet package I had to use a function called getAffy() which according to the packages' manual it :

  searches in your current or designated path file the folder with the
  GSE accession number and reads the filelist.txt file that contains the
  name of each GSM sample, in order to recognize them and join them in
  an unified AffyBatch object.

So in some way, this file is needed to create the AffyBatch object. The thing here is that there is no such a file downloaded from the getInfo() function.

Googling out the filelist.txt I reached only this result from GEOQuery manual, referring to that file. And as you can see there is a command in it eList2 <- getGEOSuppFiles("GSE11675") that is believed that it returns also the filelist.txt. I tried even with the example of the previously mentioned command, but no filelist.txt file ever fetched.

Is there any way to setup my own filelist.txt file? Does anyone know what it contains and in what form?

Finally, I found it. filelist.txt file is accessible from the FTP server of the NCBI. For example for that specific microarray is stored here

So there might be a problem at software side and they don't get downloaded.


Dealing with absence of coverage In my snakemake workflow, I use deeptools bamCoverage (in a make_bigwig rule), then computeMatrix and plotProfile (in different plot_*_profile rules) to create coverage profiles for various kinds of small RNAs. It sometimes happens that there are very few reads in a particular category of small RNAs. This results in bamCoverage exiting with an error (when I use the --skipNAs option, in order to make coverage profiles not take into account empty regions):

Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.

The fact that bamCoverage fails make things a little annoying for me, because the later the failure, the easier to handle things in my workflow: I will need to catch the error and generate a dummy output file in the make_bigwig rule, and then handle the special dummy file in the plot_*_profile rules.

I would be fine having an empty coverage file: the absence of reads (either in the bam file, or matching selection criteria in bamCoverage) is not due to an error.

Is this just not possible, by the format definition?

Did I miss an option to have bamCoverage accept making empty files?

Would this make computeMatrix fail?

I'll fix the bamCoverage issue in the 2.6.0 release (it may already be fixed since I rewrite how bigWig files are made for that release). You can follow my status on the issue page. In principle you should be able to do this, it's just a deepTools limitation.

This should work fine with computeMatrix, though keep in mind that you're just going to get a big block of NA values.

Update: This error ends up coming from pyBigWig. I've just found that bedGraphToBigWig also can't make a bigWig file with now entries, so I suspect this is sort of undefined behavior. Regardless, I'll alter pyBigWig to produce a file that it can open and use in a meaningful manner.


What are some good sources of 3rd order gene expression data I need good sources of gene expression data in the form of 3rd order tensors. Typically the commonly available datasets are in the form of a matrix, for instance, $sample \times gene$ or $gene \times time$. I need good sources for 3-dimensional data, for example $gene \times sample \times tissue$ or $gene \times sample \times time$ ... etc.

I familiar with the inflammatory bowel disease  which is quite complicated:
It affects multiple sites and depending on when it appears it is different. So you can look at the GEO for microarrays and RNA-seq data of this disease and I'm sure you'll find datasets from the same patients from several tissues or time.

In this disease there is a multi'omics project that fits your needs. You can find (and download) data of high dimensionality: there are samples for the same time, tissue (but different locations), gene, over a period of time and for several omics (not only expression data).
The only problem is that there are some mistakes on the metadata.

However, I don't think they are usually described as 3rd order tensors and I doubt you'll find enough data to train a machine learning method.


Using Broad Institute RNAseq pipeline There is a Github report for Broad's pipeline on RNA-Seq. It has a docker build file and WDL scripts.

However, I don't see how any of the WDL script is being used in the docker (Dockerfile). It looks like Broad's docker runs on the standard command-line tools, but not via the WDL language.
How do the Broad bioinformaticians take advantage of the WDL scripts?

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

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


How to use the hmmsearch on prototypic repetitive sequences Repbase Update database? I want to use the hmmsearch proposed in Convergence of retrotransposons in
oomycetes and plants by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. After the help from terdon I managed to create the correct profile I assumed that the command

hmmsearch aRNH_profile.hmm prirep.ref > output.out

would be the correct command. I downloaded the compressed database RepBase22.11.fasta.tar.gz from here in the fasta format.
However when unpacking I get the following tree

.
├── angrep.ref
├── appendix
│   ├── angapp.ref
│   ├── athapp.ref
│   ├── bctapp.ref
│   ├── cbrapp.ref
│  ....
├── athrep.ref
├── bctrep.ref
├── cbrrep.ref
...

which file is the correct fasta file? I assumed prirep.ref  might be correct file but the output shows

[No hits detected that satisfy reporting thresholds]

RepBase explain what these files are on their website. It is usually a good idea to check the documentation when downloading files so you know what you're downloading.

Based on the list of files described on the linked page, there is no prirep.ref file, so I assume you mean plnrep.ref which is "Other Plants". In any case, again according to that page, the files with plant and/or fungal data (which is what I assume you want; you never said) are:

athrep.ref  Arabidopsis thaliana
oryrep.ref  Oryza sativa
dcotrep.ref Other Dicot Plants
mcotrep.ref Other Monocot Plants
grasrep.ref Other Grass Plants
plnrep.ref  Other Plants
fngrep.ref  Fungi

There's also chlrep.ref for the green alga Chlamydomonas reinhardtii, and diarep.ref for the diatom Thalassiosira pseudonana. I don't know whether you would want to consider either here.

In any case, since these are fasta files, you cna combine them by simply cating them together:

cat athrep.ref oryrep.ref dcotrep.ref mcotrep.ref grasrep.ref \
    plnrep.ref fngrep.ref > plantsAndFungi.ref

The, run your search on the plantsAndFungi.ref file.


How to find genes whose expression vary the most in time-series? Hi I am reading this paper and trying to analyse the same datasets in my work:
dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data
In this paper, the authors are working on several time course gene expression datasets, including "Drosophila melanogaster embryogenesis" dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6186) . This dataset has the whole genome ( 14064 genes) of Drosophila melanogaster measured over 31 time-points.
In the dynGenie3 paper, we analyse a subset of these  14064 genes, as they write:
"We focused our analysis on the 1000 genes whose expression vary the most across the time series."
I want to know how do we measure the variance in expression across time series? What statistical tests might have been performed? I couldn't find any more information in the paper so I am assuming it might be something simple and frequently done.  Need some hints.

Maybe my thoughts are too simple, but wouldn't the computation of the variance per gene suffice for this? You really do not need a statistical test, since they are not talking about DEGs, just high-variant genes.


Tumor purity/contamination/admixture estimation Can anyone recommend a good tool for estimating the tumor content given a matched tumor and normal file for DNA NGS whole genome sequencing data or whole exome data?

Is it possible to estimate this without a normal sample as well?

I have previously estimated tumour purity with the EXPANDS an inferred tumour heterogeneity program which is designed to calculate the number of clonal subpopulations in matched tumour/normal samples.  The purity is essentially the size of the largest subpopulation identified in that sample - this is discussed in the programs FAQ.  In addition to a matched tumour / normal exome or genome sequencing you'll also need match somatic copy number as a seg file as input.  Other programs also exist for this sort of inferred heterogeneity analysis - some of which may also give you a measure of purity are: Absolute, ThetA, SciClone, CHAT, PyClone and Canopy. A more complete list looks to be here.

The only other thing I'd suggest is have an estimate via some orthogonal measure of the purity to try at least judge how the different methods are performing with you data, e.g. you might have clonality indications from cytogenetics/FISH or histology work, or perhaps FACS.  Picking samples known to be pure / impure via one or more of these might help you get a handle on how well various methods are performing.

As regards estimation without a normal sample the only program which looks like it might help is QuantumClone, but that requires more than one tumour sample form said patient either spatially or temporally distinct to perform the analysis.

Also If you have low coverage DNA sequencing CNAnorm looks like it could be useful with just a single sample.


How to merge transcript sequence with same name in a FASTA file. Suppose we have a fasta file like

>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCAT
>Seq1
AAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

I want to produce

>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCATAAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

Certainly this python script is not the most elegant solution but it will give you the desired result:

import sys

seqs = {}
with open(sys.argv[1], "r") as fh:
    curr = ""
    for line in fh:
        if line.startswith(">"):
            curr = line
            if curr not in seqs:
                seqs[curr] = ""
        else:
            seqs[curr] += line.strip()

for ident, seq in seqs.items():
    print("{}{}".format(ident, seq))

If you saved the script in seq_merge.py just run it with:

python seq_merge.py your_file.fasta > result.fasta


Issue with visualising cladogram/phylogenetic tree with multiple sequence alignment data in R? I would like to visualize tree with multiple sequence alignment.
My din.newick file is shown below,

(org1:0.49999997,org2:0.47298786,org3:28.37858179);

My din.fasta alignment file shown below,

>org1
ATGCTA
>org2
TAGCAT
>org3
AGTAGT

I have tried the following code,

library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta")

The above mentioned code is generating tree with the alignment matrix, but the tree root tip labels are not visible. It is either overlapped with matrix or it is not at all visible.  Kindly help me to fix this issue.

You were very close, you just need to supply the offset argument to msaplot().

So using your code block above with one extra argument in msaplot()

library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta", offset = 5)

yields

Thank you for supplying a minimal and reproducible example!


How does blastn sort output? I used blastn to search one a genome database for the sequences in a file.

blastn -out output.txt -outfmt 6 -query sequence_list -db genome_database -perc_identity 100

These sequences are all from one contig of the genome so the output was all hits in one contig. I then appended more sequences to the sequence file and used blastn again. However the weird thing is the output did not start with the same hits as the last blast. I expected if blastn iterated through sequences in a query file that it would output the original sequences first then go through the appended sequences. Since this was not the case, how does blast order its output?

I also checked if it was sorted by which contig the hit was in and this is not the case.

BLASTN should definitely report alignments in the order that the query sequences were provided.

Based on your description and my personal experience, I think the most likely explanation is some kind of mix up. It happens to all of us, especially when we're in a rush to get an answer, but even sometimes when we're being careful and disciplined.

If you can go back and recreate the result, I would be very interested to see the two pairs of input and output files. If not, I wouldn't spend too much more time worrying about it. :-)

UPDATE Alignments are reported in the order that the query sequences are provided, even in tabular mode. The order of each query's hits is sorted by bitscore or e-value, but that is a secondary sorting. I've included an example below that has two queries, seq1 and seq2 with -outfmt 6.

seq1    scaffold2       100.000 280     0       0       1       280     11131   11410   7.73e-146       518
seq1    scaffold294     88.693  283     25      4       1       279     12160   12439   6.46e-92        339
seq1    scaffold1525    85.315  286     29      9       1       277     27142   26861   3.07e-75        283
seq1    scaffold1874    83.505  291     32      10      1       280     8558    8273    1.86e-67        257
seq1    scaffold147     82.653  294     34      9       1       279     10797   10506   1.45e-63        244
seq1    scaffold2478    82.374  278     38      8       11      280     7622    7896    1.13e-59        231
seq1    scaffold3405    87.940  199     17      5       1       196     14390   14584   1.46e-58        228
seq1    scaffold395     83.060  183     23      4       106     280     12021   11839   5.40e-38        159
seq1    scaffold2853    85.965  114     13      3       1       113     20864   20753   9.16e-26        119
seq2    scaffold9       100.000 140     0       0       1       140     1121    1260    2.34e-68        259
seq2    scaffold2950    90.909  110     6       3       1       106     110     1       6.84e-34        145
seq2    scaffold3416    95.385  65      3       0       1       65      11882   11818   1.16e-21        104
seq2    scaffold3103    95.385  65      3       0       1       65      9736    9672    1.16e-21        104
seq2    scaffold5297    88.608  79      7       1       64      140     150     228     6.99e-19        95.3


finding unique headers in a fasta file using linux command line I tried to use the following command

uniq -u reference.fasta >> reference_uniq.fasta

I'd like a count of the unique headers.

If you just want the number of unique headers, you can do this:

grep '>' reference.fasta | sort | uniq | wc -l

If you want a list of the unique headers, you can do this:

grep '>' reference.fasta | sort | uniq

If you want a histogram of how many times each header occurs, you can do this:

grep '>' reference.fasta | sort | uniq -c | awk '{printf("%s\t%s\n", $1, $2)}'


Finding URL for UCSC data hub I am trying to build a UCSC track hub (a private one).

I have all my hub files on a university server in a simple folder called /hub. I am using Ubuntu 16.04 and the server is also Ubuntu (16.04.3 LTS (GNU/Linux 4.4.0-98-generic x86_64).

I have tried using hubCheck and also simply upload the hub as instructed in the UCSC instructions.

The URL I have used to indicate the location of my hub is:

http://pcxxx.xx.xx.ac.uk/home/username/projects/project-name/hub/hub.txt

However, in both cases, I get the following error:

Couldn't open http://webserver-id/home/username/projects/project-name/hub/hub.txt
TCP non-blocking connect() to webserver timed-out in select() after 10000    milliseconds - Cancelling!

I am not very knowledgeable in servers and protocols (http, ftp etc). But am I right to think that this error means that my hub is somehow not internet accessible? I'm struggling to understand what this error means. When I log onto the server, I automatically arrive in my home folder: /home/username. My hub is in /projects/project-name/hub.

Is there a way I can check if my folder location is indeed internet accessible?

Check your webserver functions from the top down.  You can do this from any WEB browser, even one that is on the webserver system itself.  For example, does this link function:
  http://webserver-id/
and provide you with a home page for that WEB service ?  If not, then there
is not apache WEB server on that system, or is mis-configured.

Next, depending upon how that Apache is configured, try this URL
  http://webserver-id/~username/
which would be a normal way to get to user directories via Apache.
This would actually be looking for a special directory in your HOME
directory, perhaps something like $HOME/public_html/
You would need to place something in public_html/ to access, for example
a tiny WEB page:

<title>my home directory</title>

hello world

If that works, then you need to place your hub files in this
public_html/ directory.


How to concatenate FASTAs with respect to chromosome? I have the following two small FASTAs, dna1.fa and dna2.fa:

# dna1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG

# dna2.fa
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

If I concatenate these two files using cat, the piped result is as follows:

$ cat dna1.fa dna2.fa > total1.fa
$ cat total1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Here, as expected, these are two files which have concatenated together. It seems this is usually what is done when one concatenates several FASTAs together in bioinformatics.

Question1: Is it possible to concatenate these "based on the chromosome"? Here is the desired result:

>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Note that these are exceptionally small FASTAs---usually one could be dealing with files of several GB in size.

Question2: The meta question here is, are there any standard bioinformatic analyses whereby the "concatenated by chromosome" FASTA would be more useful than the straightforward concatenated FASTA via cat alone?

Using the TblToFasta and FastaToTbl scripts I have posted before, combined with classic *nix utils, you can do:

$ join -a1 -a2 <(FastaToTbl dna1.fa | sort) <(FastaToTbl dna2.fa | sort) |
    sed 's/ /=/2' | TblToFasta | tr = '\n'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation

The FastaToTbl will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to join, which will join its input on lines where the 1st field is the same (by default). The -a1 -a2 means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.

Since the join will add a space between the two joined lines, we deal with that by replacing the second space on each line with a = (sed 's/ /=/2') and then pass the whole thing to TblToFasta to revert to fasta format. Finally, we use tr to replace the = with a newline so you can see where each sequence came from as you show in your example.

For longer sequences where you want to avoid sorting the input, you can do:

$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2
         }
         END{
              for(chr in seq){
                printf ">%s\n%s\n", chr, seq[chr]
              }
         }'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation

Here, the tbl format sequences (see above) are passed to an awk script which does the work:

seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2 : if the 1st field ($1, the sequence identifier) is already in the array seq, then set the value of seq[$1] to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field ($2, the sequence). If it hasn't been seen before, then set the value of seq[$1] to the 2nd field, the sequence.
END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]} : after processing everything, go through each element in the array seq and print out the data in fasta format.

However, that will mean keeping both files in memory which may be a problem for huge sequences. You can avoid that by using temp files:

$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            if(!seen[$1]++){
                print ">"$1 > $1".tmp"
            }
            print $2 >> $1".tmp"
        }' && cat *.tmp > out.fa && rm *.tmp
$ cat out.fa
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation

This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

if(!seen[$1]++){print ">"$1 > $1".tmp"} : if this is the first time we've seen this ID, print it as a fasta header into the file $1".tmp" (e.g. chr1.tmp).
print $2 >> $1".tmp" append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (rm *tmp).

As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.

Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (>chrN) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.

But I don't really know why you would want to do it. You're the one asking for it, after all! :)


Structural variant calling for low-coverage PacBio data PacBio is selling ~10x PacBio SEQUEL long reads as an upgrade to Illumina data for SV discovery.

In a clinical setting, the main requirements are proper sensitivity and specificity but also the processing of cohorts, at least families. This requires a genotyping step, such that it can be identified whether a given variant is shared by two or more individuals or whether it is not.

What are the tools of the trade for this task?

As having 50-60x PacBio reads is not an option from a economic point of view, one has to make do with 10x coverage.

There is an evaluation of PB Honey and Sniffles algorithms for low coverage PacBio datasets in this preprint and another evaluation is shown on this poster. Both reports agree that optimal is (surprisingly) combination of PB Honey and Sniffles.

Author of Sniffles have benchmarked Sniffles against PB Honey, where he shown that Sniffles performs significantly better. Take a look on this presentation (slide 15).

Another option is SMRT-SV, but I am not aware of any benchmarking.


Estimate the length of poly-A tails from randomly-primed RNAseq data So a poly-A tail is a long chain of adenine nucleotides that is added to a messenger RNA (mRNA) molecule during RNA processing to increase the stability of the molecule.

For my project, I would like to estimate the length of poly-A tails from randomly-primed RNAseq data.

Any ideas how this could be done? Are there algorithms/tools which currently do this?

Perhaps I would need to model this somehow? All suggestions appreciated

Can’t be done. If you already sequenced then I’m afraid the money is wasted (unless of course the data is good for something else).

The standard Illumina basecaller doesn’t deal well with homopolymers and will consequently call incorrect poly(A) lengths.

Amongst others, Narry Kim’s lab developed a dedicated method — TAIL-seq — to solve this problem. It’s a protocol for library prep using pulldown of biotin tagged tails (rather than random priming), plus dedicated basecalling software that uses spike-ins for calibration.

Be warned though: the protocol is reportedly technically challenging, and in order to use it you need an Illumina machine that gives you access to the raw image files before basecalling. Up-to-date Illumina machines won’t give you these files, full stop. So you’ll need to find a machine that runs an outdated firmware.


Extracting some part of a list I have two normalized gene expression values (log2 of cpm)

> head(biomarker[,1:5])
               A2        A3        A4        A6        A7
A2M    12.4071618 12.631601 12.889748 11.842511 11.574134
ABCB11  0.5151957  6.176788  5.414905  7.134915  1.247590
ABCC2   6.7244303  6.487794 10.178615  6.384132  6.089462
ABCC6   6.3977122  5.823627  7.370761  7.397091  6.071587
ABCF1   9.1609847  9.893258 10.116638 10.225520 10.838486
ABCG2   7.6841874  7.738293  6.833152  6.888041  6.030032
>

And

> head(immune[,1:5])
              A2        A3        A4        A6        A7
A2M    10.904748 11.388404  9.910614 11.513439 12.963609
ABCB11  5.011380  6.359443  7.145992  8.451947  7.722605
ABCC2   5.040461  6.477014  3.873996  6.409777  9.133971
ABCC6   7.798441  7.601848  9.948072 11.628533 12.701460
ABCF1   8.553597  8.615425 11.145903 10.289098 11.444140
ABCG2   6.224294  5.629375  8.293416  7.979859  8.603793
>

> dim(biomarker)
[1] 719  56
> dim(immune)
[1] 719  56
>

I have matched samples in these experiments and same gene too. I could roughly say in both experiments everything has been constant. However, I want to do Student's t-Test to find genes that they had inconsistent gene expression patterns between the two data sets.

By this function I am getting a list but I don't know how to extract the names of genes with p-value > 0.05 or whatever cut-off

> f <- function(x,y){
+     test <- t.test(x,y, paired=TRUE)
+     out <- data.frame(stat = test$statistic,
+                       df   = test$parameter,
+                       pval = test$p.value,
+                       conl = test$conf.int[1],
+                       conh = test$conf.int[2]
+     )
+     return(out)
+ }
> t=sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Any help please?

  str(t)

    [list output truncated]
     - attr(*, "dim")= int [1:2] 5 719
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:5] "stat" "df" "pval" "conl" ...
  ..$ : NULL
    >

> sapply(t, "[[", "pval")
Error in FUN(X[[i]], ...) : subscript out of bounds
>

> pvalueGenes <- sapply(colnames(biomarker),function(i) t.test(biomarker[i, ],immune[i, ], paired = TRUE)$p.value)
Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mean.default(x) : argument is not numeric or logical: returning NA
2: In if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :

 Show Traceback

 Rerun with Debug
 Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed

From your example we are comparing samples not genes, I think we need to run the test row wise, something like this:

pvalueGenes <- sapply(colnames(biomarker),
                      function(i) t.test(biomarker[i, ],
                                         immune[i, ], paired = TRUE)$p.value)

names(pvalueGenes)[ pvalueGenes < 0.05 ]


Find CRISPR PAM-sites with python I am a college student and just started bioinformatics. I am trying to write a script (not for school) that finds all potential NGG and TTTN protospacer adjacent motif sequences in a genome string. I wrote the following script:
TargetFile = open('file.fasta')
readTF = TargetFile.read()
useTF = "".join(readTF[66:].split())

#no_of_basepairs = len(useTF)

sites = {"ngg_strand" : [],
         "ngg_antistrand" : [],     #location counted from reverse direction
         "tttn_strand" : [],
         "tttn_antistrand" : []}    #location counted from reverse direction

for i in range(0, len(useTF)-2):
    gg_mer = useTF[i:i + 2]
    ttt_mer = useTF[i:i + 3]

    if gg_mer == "GG":
        sites["ngg_strand"].append(i)
    elif gg_mer == "CC":
        #look for reverse complement, record location on antistrand from reverse direction
        sites["ngg_antistrand"].append(len(useTF)-i-2)
    elif ttt_mer == "TTT":
        sites["tttn_strand"].append(i+1)
    elif ttt_mer == "AAA":
        # look for reverse complement, record location on antistrand from reverse direction
        sites["tttn_antistrand"].append(len(useTF)-i-2)

print(len(sites["ngg_strand"]))
print(useTF.count("GG"))

At the bottom I used the .count() method to compare results. I noticed that the method returns a smaller number than my code. I believe this is due to the .count() method not counting overlap, whereas my code does. I.e. in the string "NNNGGGGNNN" my code would count 3 occurences of the "GG" substring, whereas the .count() method returns 2 occurences.
My Questions:

If I want to find all possible PAM sites, which way of counting is correct? I believe my way would be correct, since as far as I know the RNA molecule doesn't advance in "triplet-steps". However, I am very much new to this and might be missing something.

If the reverse complement is found, is the location usually recorded from the same direction as on the strand, or in the reverse direction (as I implemented in the code)?

Like I said, I am new to this so if you see any other glaring problems I would appreciate any and all help. :-)

Thank you!

The easiest solution is finditer
use re

# Compile you search term where [A|T|G|C] means A or G or C or T
p = re.compile("[A|T|G|C]GG")

iterator = p.finditer(useTF)
for match in iterator:
    print(match.span())

** Example Output**
(0, 2)
(22, 24)
(29, 31)

Easy does it. Its worth looking at the documentation for re here, https://docs.python.org/3/howto/regex.html for alternative approaches using match or search below. finditer is definitely the solution in context.

I support your coding approach because it provides lots of flexibility. For example, you might want to target GG or TTT around a given locus and you can easily modify the code to accomodate this. Someone elses code would be less flexible.
In summary, I kinda like your substrings and do like the count method. I'm not a fan of your == in this instance and your C-style matching should be substituted for reg-ex which in Python is implemented via re.
The basic issues are,

Biopython would be used via AlignIO or SeqIO. This will allow any given format to be analysed. Yours approach is okay but prone to bugs (below). [probably a minor issue]
Your sequence is in some sort of weird format, like extended phylip format and it is not fasta (that wouldn't work). It needs verifying because removing the first 65 characters might remove bits of sequence [probably a minor issue]
Chopping of -2 from the length ... I dunno I'd check out -1 [minor]
C-style matching. I don't like this. I would definitely use re. One solution uses match or search, but the solution above uses finditer, Reg-ex will allow to search for NGG and TTTN, which you are not doing
I think you are looping a locus of 3 and 4 bp and I'm not sure how the C-style == would match this because thats 2 and 3 bp respectively. Just thoughts. You will have 'edge effects' at the beginning and end of sequences, e.g. TTT 3' will match but that isn't TTTN.
I like your count method, but I think finditer is better because it gives you location information
The dictionary is sort of fine. Personally, I would dump the dictionary into a pandas dataframe This way the results are output as e.g. CSV, spreadsheet flat-files for your colleagues to view.
the if and elif would be placed in a subroutine, this is good practice to convert to an object orientated approach.


Tool for calibrating base quality scores? A given sequencing machine assigns a 'certainty' to each base call in the form of a per-base quality score (FASTQ format).

I have reads from several machines aligned to the corresponding reference (BAM format or similar). Ignoring SNPs and misalignments (for simplicity), mismatches should represent sequencing errors.

I'd like to parse through the FASTQ and BAM file and evaluate how closely the per-base quality score reflects 'true' sequencing errors.

Is there an existing efficient, user-friendly, open source, well documented tool to do this? Failing that, is there a paper doing this that claims to have made the tool freely available?

Another method would be to look at bases corrected using kmer analysis or similar. Again, the per-base quality score can be compared to the 'true' sequencing errors.

Any tool to do that? I guess it could be baked into an existing error correction tool.

Do short read aligners actually make use of per-base quality score?

Many thanks,

I think what you are looking for is GATK's Base Quality Score Recalibration.  Their documentation is great so I won't rewrite it here.

  BQSR stands for Base Quality Score Recalibration. In a nutshell, it is
  a data pre-processing step that detects systematic errors made by the
  sequencing machine when it estimates the accuracy of each base call.


Filter out outliers of the scRNA-seq (heterogenous cells) I am new to data science. I have a dataset of single-cell gene expression from multiple cell types in C. Elegans. The dataset is from the paper Comprehensive single-cell transcriptional profiling of a multicellular organism

My main question is, Which approaches should I use for filtering out bad cells in this case when we have multiple cell types in the dataset?

So far I tried to filter out genes that have too high mitochondrial genes content following the Bioconductor “simpleSingleCell” workflow.

However, the tutorial specifically says that the method of filtering out based on mitochondrial genes most probably will not work when the dataset has multiple cell types:

  Analyzing all cell types together would unnecessarily inflate the MAD and
      compromise the removal of low-quality cells, at best; or lead to the
      entire loss of one cell type, at worst.

Any suggestions would be greatly appreciated.

From what I known, there is no clear consensus in the field and it depends on the type of cells you are interrogating.

However, if mitoRNA/endogenousRNA ratio does not fit for your purposes, other option is to check at the total number of genes/transcripts detected in each cell. In this way, you can filter out cells that have considerable less genes/transcripts detected than the rest since this may be indicative of bad quality cells for any reason (apoptosis, RNA degradation, sequencing itself, etc). For instance, to put a threshold of >1000 genes/transcripts in a cell to be considered for further analysis.

Consider also to take a look at the number of mapped reads for each cell, as cells with low mapped reads respect to the rest could be potentially problematic.


Which module to select for Pathway analysis based on Module trait correlation and pvalue? I have a total of 35 tumor samples classified into 4 subtypes. Subtype A, B, C, and D. I have RNAseq data. I'm interested in identifying modules related to each subtype with co-expression network analysis. I have read the WGCNA tutorial and create a trait.csv file like below:

Showing here some rows of the trait.csv file.

SamInfo <- read.csv("trait.csv")

Samples = rownames(datExpr);
traitRows = match(Samples, SamInfo$Samples);
datTraits = SamInfo[traitRows, -1];
rownames(datTraits) = SamInfo[traitRows, 1];

nGenes = ncol(datExpr);
#15648L

nSamples = nrow(datExpr);
#35L

#calculate eigengenes (1st principal component) of modules
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

#Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")",
                   sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));

#Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits),
               yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE,
               colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE,
               cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships"))

The Correlation between Subtypes and Modules Heatmap looks like below:

For functions identification (Pathway Analysis), from the above output for each subtype, should I select only the modules that were significantly positively correlated with subtypes, or Can I also use the significantly negatively correlated modules?

Is it good to use the modules that were not significant for pathway analysis?

In case if I can use positive and negative correlated modules, Can I combine those modules and use those genes for Pathway analysis?

If you can run pathway analysis programmatically on all modules, yes, do it on all modules. Otherwise, select the modules with strongest associations with each subtype, both positive and negative. I would keep the modules separate because they may have very different biological functions.


RNA-Seq: clustering/treatment of genes with low expression I have some RNA-Seq data from leukaemia patients. I want to do unsupervised clustering on them with some other published leukaemia RNA-Seq data and see how they cluster. There are a few problems I encountered while doing this.

I read mix messages of whether using log2(TPM) or rlog(counts) (e.g. Deseq2 rlog or limma voom transformation) for clustering. Which one should I choose?
If I am filtering out genes with low counts, should I do it prior normalisation to library size or after?
I tried using filterByExpr from edgeR for filtering but it removed many of the genes that are only expressed by subtypes with small sample number. If I am to filter the counts "manually", is there a recommendation of how the threshold should be set?
If I should use log2(TPM) for unsupervised clustering, can I treat low TPM (e.g. <2) as 0? As this paper suggested TPM < 2 basically means no expression, is it alright to treat them as 0?

rlog(normalized counts) is going to be more robust than log2(TPM), so use it instead.
Do it afterward, keeping the low counts can be helpful for library-size normalization.
This is very much more of an art than a science, though have a look at the kOverA() function from the genefilter package. That will give you a bit more fine-grained control on filtering.
You'd be better off using something like zFPKMs if you want to filter by "expressed or not", since there's no really great definition of that.


different results coming from biomart online and biomaRt R library I am playing with biomart in order to get the list of genes associated to a specific gene ontology (GO) term.

So, if I use the webservice as in this query:

I get at least 50 genes. If I download the file I can gen an exact count of them:

$ cat mart_export.txt | wc -l
347

Now, I can do the same operation by using the R library biomaRt. Specifically:

> ensembl <- useEnsembl(biomart="ENSEMBL_MART_ENSEMBL", dataset='hsapiens_gene_ensembl')
> gene <- getBM(attributes = c('external_gene_name'), filters = 'go', values = 'GO:0030098', mart = ensembl)
> unique(gene)
  external_gene_name
1               SPI1
2              PRKDC
3                TOX
4               LY6D
5               CBFB
6              IKZF1
7               RELB
8             FLT3LG

By doing so, these are all the genes I retrieve. The two results are totally different and I don't know why it. All these eight genes are found in the biomart file (I grepped it, to be sure).

I believe this is not an issue related to databases versions since I checked:

>listMarts()
               biomart               version
1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 94
2   ENSEMBL_MART_MOUSE      Mouse strains 94
3     ENSEMBL_MART_SNP  Ensembl Variation 94
4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 94

and:

>listDatasets(ensembl)
...
38 Human genes (GRCh38.p12)
...

Both are the same version used by biomart online.

What is going on? Is it possible that biomart online is considering GO relations that are not specified in the R version?

Sorry, this was my mistake in the last question. To search down the ontology, rather than just for the specific association with a term, biomaRt needs a different filter: 'go_parent_term'.

Try:

gene <- getBM(attributes = c('external_gene_name'), filters = 'go_parent_term', values = 'GO:0030098', mart = ensembl)


Compute the significance of the overlap between 2 or more gene sets I was able to compute the significance of the overlap between 2 gene sets using the cdf function of scipy hypergeometric distribution.
I wish to be able to perform the same calculation for more than 2 gene sets; should I use the multivariate hypergeometric distribution cdf function for that?
Are there any websites that provides the same calculations over gene sets so I can validate my results?

This is a shot at it, first an example dataset:
import matplotlib.pyplot as plt
import numpy as np
import functools
from matplotlib_venn import venn3

# define universe
uni = ["gene"+str(i) for i in range(1000)]
# some overlap
gs1 = uni[250:300] + uni[900:950]
gs2 = uni[:300]
gs3 = uni[250:500]

The ever amazing venn diagram:
venn3([set(gs1),set(gs2),set(gs3)],set_labels=["gs1","gs2","gs3"])

Then a function to draw a set with length equivalent of each set, randomly from the universe and find length of intersection (all 3):
def sim_intersect(uni,set_lengths):
    randomsets = [np.random.choice(uni,n) for n in set_lengths]
    return len(functools.reduce(np.intersect1d,randomsets))

We run this 1000 times:
permuted_values = [sim_intersect(uni,[len(gs1),len(gs2),len(gs3)]) for i in range(1000)]

plt.hist(permuted_values,bins=range(50))

The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this post:
(sum(np.array(permuted_values)>obs_n)+1)/(1000+1)
0.000999000999000999


what is the concept of N-channel geneChip? I can understand two-color and one-color(Channel) gene chip; but there is another type of gene-Chip called "N-channel chips" ; I do not know how can several (> 2) samples be hybridize in one chip. I cannot find some reference for the a material for this concept. See the following code that it has two channels "E“ and "Eb"; what does that channel mean? The expression matrices seem to have the same dimension. In my understanding, the two-color or two-channel microarray is to combine two samples (treated vs. control) in a chip, and the log(ratio) is used to express the differential expression.

AE1548 = getAE("E-MTAB-1548",
               path="/Users/zhang/Documents/2020/GEOsepsis/Data",
               local = TRUE, sourcedir = "/Users/zhang/Documents/2020/GEOsepsis/Data",
               type = 'full')
AE1548raw = ae2bioc(mageFiles = AE1548)#established from raw data; multi-channel data
AE1548expr1 = assayDataElement(AE1548raw,"E")
AE1548expr2 = assayDataElement(AE1548raw,"Eb")
  AE1548expr1[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1    273812.0     97733.5      205330
2        41.0        26.5          50
3        42.5        25.0          48
4        40.0        23.5          49
5        38.5        24.0          53
  AE1548expr2[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1          25        18.0          36
2          25        18.0          37
3          25        18.0          38
4          25        19.0          37
5          24        18.5          37

This is a mere long comment, which may help, but I do not know the details required for the answer.

Several (> 2) samples hybridized in one chip

A traditional microarray has two parts, the mobile target, which is the labelled sample, and fixed to the filter paper, glass or doped silicon is the probe, which is not labelled and is in vast excess. You can label the target fluorescently with one of a variety of fluorophores, or with different radioisotopes, which have different energies. Nobody uses radioisotope microscopy for this except for the most obscure method, say nuclear run off assays, but I mention this just to cover all bases.

Whereas you can buy whatever fluorophore you wish, the detecting machine needs to be set up for it with lasers, beam splitters and filters. Most florescent confocal microscopes use a LED lamp, a splitter and a loads of cubes of filters. In the case of flow cytometry, you can get over 16 fluorescent channels at the same time. However, deconvolution becomes necessary when one channel "bleeds" into an other —I am not sure how one deals with it with a microarray, I suspect before the step to remove regional biases (the step with LOESS regression or similar).

Affymetrix Gene chip (R)

Now, in the case of an Affymetrix Gene chip, I do not know if there is a multichannel machine, but it could be custom build. The spots on an Affymetrix chip are 3-5 µm, so regular microscopy optics still work. Although the silicon base is not transparent like a home made slide or an Agilent microarray (if these are still made). Hopefully, someone else can answer about Affymetrix machinery!

Statistics

I am not familiar with why or how you would analyse a 3+ channel experiment. One guess could be to do a dose-response experiment. But there will be some problems:

The dye flip correction would be an issue for greater than 2 samples
Channel deconvolution may be required for dyes with overlapping spectra.


At what stage of a transcriptome assembly is it better to perform read contaminant filter? I'm trying to assemble a bivalve transcriptome. Since bivalves are filter feeders, their transcriptomes tend to be highly contaminated by bacteria, algae and whatnot. Since I pooled several transcriptomes, I have a high amount of reads (>2B reads). I also have a not-so-high quality reference genome.

Now, I'm not sure when to perform a BLAST-based filter on these reads:

Do I filter reads with BLAST before assembling them?
Do I assemble the transcriptome and then filter out contaminant contigs?
Since I have a reference genome, do I simply not care about contaminant reads because they will not align to it?

Having this amount of reads, I'd be inclined to first filter out those reads that don't have a match against Metazoa, and then to perform a second filter after assembly in order to retain only those contigs that match against Lophotrochozoa.

If you know resources on how to manage read contamination, feel free to share them.

First, to answer your question about mapping to a low-quality reference:

1. Mapping

For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome assembly does not contain any scaffolds from contaminating species like bacteria, algae, et cetera.

2. Reasons to not BLAST your raw reads

There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in #4, then reassembling, may work better.

3. Using multiple individuals of highly heterozygous species in assembly.

Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.

4. A strategy to filter transcripts based on mapping and BLAST

Given #1, #2, and #3 above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.

I would try this:

Assemble the transcriptome of each individual separately.

This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.

Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.

Use somthing like cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta.

Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.

Use samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam to only keep the transcripts that did not mapped to the reference. -f 4 is the magic flag here that keeps things that didn't map.

Make a list of all the transcripts that did not map to the reference...

samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt

...and get only the transcripts from your fasta file that are potential contaminants.

seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta

BLAST all of the unmapped.fasta and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this known_contaminants.fasta
Map all reads from all individuals to the known_contaminants.fasta.
Get all the reads that map to known contaminating transcripts, known_contaminants.fasta using:

samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt

Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
Reassemble the filtered fastq files.

This strategy has the following benefits:

Does not require weeks/months of compute time to BLAST 2 billion reads.
Ensures that all reads from the reference genome are left untouched.
Ensures that all reads from known contaminants are removed.
Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
Leverages contaminants from different individuals to filter reads from all individuals.


Generating 3D coordinates error I have been trying to convert a multi entry sdf file to pdb using openbabel. Everytime I check the Generate 3D coordinates box I end up with .ENT files instead of .PDB files. What am I doing wrong.
Thank you.

An ent file is a pdb file. Just change the extension. PDBe for example provides their PDB files as such.


Moving file based on their names I have a list of vcf files; I also have a list of names in a txt file like

LP6005409-DNA_F01
LP2000325-DNA_A01
LP6005409-DNA_E02
LP6005500-DNA_C03

What I have in this txt files are part of name of each vcf files

I want to move .vcf files having these in their name to a separate folder

I done

[fi1d18@cyan01 TRG]$ for i in *.vcf
> do if grep -q $i 1.txt;
> cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
> done
-bash: syntax error near unexpected token `done'
[fi1d18@cyan01 TRG]$

How I can complete this task?

Thank you

First of all you are getting an error because you are missing a then in the third line of your code and then a fi to close the if. It should be:

for i in *.vcf
do if grep -q $i 1.txt; then
    cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
fi
done

And I guess the one-liner below achieves what you are trying to do:

while IFS= read -r line; do mv "$line"*.vcf some-directory; done < file-list.txt


Converting aligned fasta to plink ped/bed I have an alignment of multiple sequences in a FASTA file (output from MAFFT), for which I would like to simulate a phenotype using plink, but for that I need to have my alignment in a PED file or plink BED. There's plenty online about converting in the other direction, but the only mention I found of converting FASTA to BED is about aligning a raw FASTQ file for use in plink - my data is already aligned! How can I do this conversion?

Note that How to convert FASTA to BED is about the UCSC BED format, not the plink binary PED.

An alignment can be the result of two slightly different analyses:
There is multiple sequence alignment (which is what you get from MAFFT) where sequences are aligned so that similar regions are on top of each other. This may require introducing indels (insertions / deletions) if a particular region is absent in some of the sequences.
Then there is alignment to a reference sequence. This is usually how you get variant data and most of the tools deal with this use case. Therefore, in the link you provided the solution includes aligning the data to a reference sequence. Usually the refernce sequence is much larger (e.g. the human genome) and after alignment the result is a file (.bam) that tells you where the query sequences match in the reference.
I did a quick search for ways of converting a multiple sequence alignment (MSA) to VCF. There is a tool called msa2vcf in the Jvarkit collection of utilities that can do this. The example is for the CLUSTAW format but FASTA is accepted as well.
You shoud check that indels were treated correctly because these are most likely to cause trouble.
Then you can simply convert the VCF to PLINK format  using PLINK.


What are the p01sum and p0prop parameters in godon? When using branch-site model in Godon, model optimization results include the following parameters: p01sum and p0prop. What are those and how do they relate to p0 and p1 in the branch-site model?

Disclaimer: I'm the author of godon. The question above comes from the real communication with a software user.

In the godon implementation of the branch-site model I performed the following reparametrization:

$p_{01sum} = p_0 + p_1$

$p_{0prop} = \frac{p_0}{p_0 + p_1}$.

This has the following advantages:

Both parameters have a well defined range $(0, 1)$.
There is less dependency between the two. E.g., $p_2=1-p_0-p_1$ (proportion of sites under positivie selection) does not depend on $p_{0prop}$.
Both properties are very helpful for both likelihood maximization and MCMC.

In case you are interested in $p_0$ and $p_1$ there is a straightforward way to go back to the original parameters:

$p_0=p_{0prop}*p_{01sum}$

$p_1=p_{01sum}-p_0$.

P.S. I also added this information to the tutorial.


Seurat VlnPlot presenting expression of multiple genes in a single cluster Seurat VlnPlots are most commonly used to visualize differences in any given gene expression across multiple clusters or cell types. For example:

VlnPlot(object = MouseCellAtlas, idents = c("T cell","Neutrophil","Erythroblast","Monocyte","Macrophage"), features = 'Fzd9')

But suppose I want to view the expression of several genes in just idents = "T cell". I can do this to generate three separate plots:

VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

But ideally I'd like this second plot to look just like the first plot, except with only a single tissue on the legend ("T cell") and multiple genes in the various colored violins. Unfortunately, I am not aware of any simple way to combine these violins into a single plot area.

This doesn't do what I need:

VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

Is there built-in functionality in Seurat for generating a VlnPlot for multiple genes in a single tissue in a single plot? Or can anyone suggest a workaround?

To do so one workaround it to have your data in "long format" and then use the column that holds the "gene names" as the x variable while plotting.

You can use FetchData() to extract data from a Seurat object. VlnPlot's default is the data slot (of the active assay if using Seurat v3 I suppose). And you can specify which cells and genes to retrieve.

selected_cells <- names(panc8$celltype[panc8$celltype == "gamma"])
data <- FetchData(panc8,
                    vars = c('FZD9','CTNNB1','APC'),
                    cells = selected_cells ,
                    slot = "data")
> head(data)
           FZD9   CTNNB1      APC
D101_5        0 0.000000 0.000000
D101_43       0 2.007853 1.001958
D101_93       0 0.000000 0.000000

Melt() would transform your data into "long format". By not specifying any arguments, we push for all of the info in the three variables to be gathered in two columns:

long_data <- melt(data)
No id variables; using all as measure variables
> head(long_data)
  variable value
1     FZD9     0
2     FZD9     0
3     FZD9     0

> tail(long_data)
     variable value
1870      APC     0
1871      APC     0
1872      APC     0

ggplot2 is used to add "violin" and "jitter" layers. You can customize the output to look (exactly) like the VlnPlot() output.

ggplot(long_data,
       aes(x = variable, y = value)) +
  geom_violin() +
  geom_jitter(size = 0.1)

And here is the same graph generated by VlnPlot():

There are no "violins" as the counts are almost entirely zeros. And see how different the ranges are in the y axes of the VlnPlot. So the "tweak" I have presented here would only work for genes that are expressed at similar levels / similar ranges.


Targeted NGS, up to 99% of reads have been marked as duplicates Currently I'm performing whole analysis (pipeline from *.fastq to *.vcf) of 41  samples (targeted NGS). I rely on GATK best practices, however with some modifications. I decided to use the following tools:
#mapping
bwa mem (with mem - alternate haplotypes)

#preprocessing
RevertSam (Picard)
MergeBamAlignment (Picard)
MarkDuplicates (Picard)
BaseRecalibrator (GATK)
ApplyBQSR (GATK)
AnalyzeCovariates (GATK)
SortSam (Picard)

#call, filter and annotate variants
HaplotypeCaller (GATK)
GenotypeGVCFs (GATK)
snpEff, SnpSift

Nevertheless, after marking the duplicates with MarkDuplicates about 99% of reads have been marked as duplicated, hence the coverage decreased from several hundreds to literally up to five. So I decided to omit this step.
Have you ever encountered a similar problem? What is your opinion about this?
Best Regards,
Adam
EDIT
The Agilent Haloplex custom target enrichment kit (HaloPlex Target Enrichment for Illumina Kit) has been used to prepare the libraries, sequencing done using Illumina HiSeq 2500. Targeted NGS of all exons (only) of 90 genes.
Here I provide how the final bam look like. This is the result without duplicate marking, otherwise there're only several reads per exon...

If you're doing target capture, especially of a small region, it is entirely possible that there are only five unique reads. Target capture is molecularly finicky and sometimes the entire region is not captured.


Consensus sequence from SAM or BAM file? I am trying to perform reference-based assembly. Most of the tutorials teach how to create a bam file and view alignemnts in IGV or Tablet. But, I want a assembled genome sequence in fasta format. How can I do that ? Is there any tool or linux syntax available for this purpose ?

A few years ago I used AlignGraph (paper and code) for a task like this. I gave AlignGraph my raw reads, so I'm not sure if it accepts SAM/BAM input or if you have to go back to the original Fastq reads.

I'm not sure if any related/improved tools have been published since AlignGraph. Internet searches for "reference-guided assembly" return mostly transcriptome assembly tools. It might be worth looking for articles that cite the AlignGraph paper.


Show presence of known mutation in RNA-seq data We have RNA-seq fastq data from control (WT) patients and a patient with a point mutation at a known location in one gene.
I'd like to retrieve the reads aligning to that gene and show the presence of the mutation.

I can think of 2 strategies:

A)

Alignment and visualization, based on GATK:

Run STAR 2-pass
Extract alignments of interest, e.g.:

samtools view input.bam "Chr10:18000-45500" > output.sam
Visualize:

IGV
IGB
ASCIIGenome (see also this SE answer)

How can I count the reads showing mutation vs wild type?

B)

The second approach would simply look for the count numbers mapping (e.g. salmon) to the transcript ID(s) corresponding to mutation and WT, based on this SE answer.

Is there a better way to do this?


Creating a new reference genome B by changing genome A with mismatches from BAM file without long reads This question is slightly related to this one:
Improve a reference genome with sequencing data

Only in my case, I have a starting genome reference A, of a short phage genome, and reads from a related genome B, for which I haven't been able to find a reference.

I took the Illumina 2x75bp paired end reads from genome B and ran Velvet on it (I ran Ray, but it didn't complete). It gives me a fragmented collection of contigs. The contigs align 99-98% to reference genome A. I don't have long reads.

My question is: with the short reads and the contigs that I have, can I transform the reference genome A into a reference genome B?

PS: no, I don't have long reads, and there is no budget for them. The budget for long reads is zero and will remain zero. Even if we could sequence long reads for free, we would still not do it because of the lab work.

Tools that I have tried so far:

ococo:

$HOME/ococo/ococo -i contigs.bwamem.bam -c 1 -f $HOME/Genomes/T2_phage/starting.fasta -F consensus.fasta --vcf-cons - -V variants.vcf

It produced a consensus.fasta, and checking how much it changes from the reference, it changes. There are one or more changes in about 33% of the 100bp windows:

diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep '\ \|\ ' | wc -l
1272
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep -v '\ \|\ ' | wc -l
419

You should be able to use something like Pilon to convert from one genome to another, assuming the changes are all fairly small local changes (but that's almost never the case).

I find your resistance to long-read data interesting. It's important to consider the cost of bioinformatics as well. Sequencing a phage genome with long reads will take less than an hour, and probably produce a few reads that cover the entire genome. With existing short-read data, those full-length reads should be able to be fixed up within a day.

I spent about half a year trying to fix a eukaryotic genome assembly that was assembled only from short reads (using Pilon, and L-RNA-scaffolder), and I was only able to reduce the number of contigs down by about half (and even at that, I wasn't convinced that contigs were being properly joined). With long reads and another couple of weeks of computer work, the contigs dropped down to about a tenth of their original number.

If you can use a pipette and have about half an hour spare, you can do rapid long-read sequencing on the MinION.


Batch alignment of inconsistently named Fastq files I have numerous gzip-compressed paired-end Fastq files with ChIP-seq data that I would like to align with bowtie2. I confirmed that bowtie2 takes .gz files, although I do see this warning when I run for two samples:

  Warning: gzbuffer added in zlib v1.2.3.5. Unable to change buffer size
  from default of 8192.

My files are named like so.

10_S10_L002_R1_001.fastq.gz  12_S12_L002_R1_001.fastq.gz  2_S2_L001_R1_001.fastq.gz  4_S4_L001_R1_001.fastq-003.gz  6_S6_L001_R1_001.fastq-002.gz  8_S8_L002_R1_001.fastq.gz
10_S10_L002_R2_001.fastq.gz  12_S12_L002_R2_001.fastq.gz  2_S2_L001_R2_001.fastq.gz  4_S4_L001_R2_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz  8_S8_L002_R2_001.fastq.gz
11_S11_L002_R1_001.fastq.gz  1_S1_L001_R1_001.fastq.gz    3_S3_L001_R1_001.fastq.gz  5_S5_L001_R1_001.fastq-002.gz  7_S7_L002_R1_001.fastq.gz      9_S9_L002_R1_001.fastq.gz
11_S11_L002_R2_001.fastq.gz  1_S1_L001_R2_001.fastq.gz    3_S3_L001_R2_001.fastq.gz  5_S5_L001_R2_001.fastq-003.gz  7_S7_L002_R2_001.fastq.gz      9_S9_L002_R2_001.fastq.gz

The bash commands I tried are below. When I run this test command:

for i in $(ls *.fastq*.gz | grep -v 6 | rev | cut -c 16- | rev | uniq); do
    echo $i;
done

...I get this output.

12_S12_L002_R1_001.fastq.gz  12_S12_L002_R2_001.fastq.gz

When I run this:

for i in $(ls *.fastq*.gz | rev | cut -c 16- | rev | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}R1_001.fastq.gz -2 ${i}R2_001.fastq.gz | samtools view -bS - >${i%}.bam
done

I do see one BAM for my paired end files. But as I see the second set of files which are named:

6_S6_L001_R1_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz

if I run the same script over them, I would get an error as the loop can't read the paired end file.

So the tedious way is I can make a separate folder for these sets of data with fastq-002.gz and fastq-003.gz and can run the script by tweaking the cut parameter.

But I would like to know what is the programmatic way of resolving the issue.

A very simple approach, but one which will work with the names you show, is to remove any occurrence of 00N (where N is any number) from right before the .gz. To do so, edit your script to:

for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do
  bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \
    -1 ${i}R1_001.fastq.gz \
    -2 ${i}R2_001.fastq.gz |
        samtools view -bS - >${i%}.bam"
done


UCSC packages/tools ucsc-netchainsubset and ucsc-fatotwobit not working on GNU/Linux I have installed some ucsc tools on my Mac computer and it worked fine:

Installation

conda create --name ucsc-netchainsubset377 --channel bioconda --yes ucsc-netchainsubset=377

Conda environment activation

source activate ucsc-netchainsubset377

Checking the tool's help

    netchainsubset
netChainSubset - Create chain file with subset of chains that appear in the net
usage:
   netChainSubset in.net in.chain out.chain
options:
   -gapOut=gap.tab - Output gap sizes to file
   -type=XXX - Restrict output to particular type in net file
   -splitOnInsert - Split chain when get an insertion of another chain
   -wholeChains - Write entire chain references by net, don't split
    when a high-level net is encoundered.  This is useful when nets
    have been filtered.
   -skipMissing - skip chains that are not found instead of generating
    an error.  Useful if chains have been filtered.

The same commands on GNU/Linux produced the following error:

netchainsubset
netchainsubset: command not found

even though the tool seems to be installed successfully:

Downloading and Extracting Packages
ucsc-netchainsubset- | 310 KB    | ############################################################################ | 100%
mysql-connector-c-6. | 4.4 MB    | ############################################################################ | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

The same happened with ucsc-fatotwobit and ucsc-liftup and I would not be surprised if this behaviour affected the whole suite of ucsc programmes installed through conda.

Anyone can suggest a solution?

Update 5 minutes after posting

Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.

ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt

Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.

Cheers terdon.

Cheers.

Short answer: respect the case of the tool name

Longer story below

Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.

ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt

Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.


Detecting portions of human proteins with high degree of microbial similarity I'm a newcomer to the world of bioinformatics, and in need of help solving a problem.

My goal is to take a list of human proteins, and identify segments (13-17aa in length) with a high degree of similarity to microbial sequences. Ideally, I would like to start with list of FASTA sequences, and have an easy way to generate an output of the corresponding high similarity segments of each protein.

Are there existing tools or software that I should be aware of that will make my life easier?

Thanks in advance.

Sounds like precisely the job BLAST was developed for. Now, which flavor will depend on what you want to do and what data you have available. Some options:

PSI-BLAST: this is usually the best choice if you are trying to find protein homologs. It works by building a hidden markov model describing your query sequence and using that model to query a database of proteins. The advantage is that it is run in multiple iterations, giving you the chance to add or remove results (so you add the ones that are true positives and remove false ones), eventually building a pretty good model of your protein. This is far moer powerful than a simple homology-based approach since proteins work via protein domains and simple homology is not as important as specific conserved functional residues.

For this, go to the NCBI protein blast page and select PSI-BLAST:

BLASTp: Simple protein-protein blast. It will identify homologous proteins based on sequence similarity. Whether or not that also implies functional homology is not that simple and will depend on each case you investigate.

As above, go to the NCBI protein blast page, but this time use the defaults.
tBLASTn: this is a tool that takes protein sequences as input and compares them to a database of DNA which is dynamically translated in all 6 possible reading frames. Very good for finding homologous sequences when you don't have well annotated protein information for the target species. It has the benefit of being more sensitive and able to find more distant homologies than basic nucleotide BLASTn and the go-to approach when your target species is distant and not well annotated.

NCBI's tBLASTn page.

All of these can be run online through the NCBI's BLAST page. If you want to investigate hundreds of proteins, I suggest you install blast locally. You can then either download the relevant target sequences from NCBI and rebuild the blast database locally (if so, I suggest you ask a new question about how to do that) or, use NCBI's remote blast client which lets you use a locally stored query file and will run blast on the NCBI's servers.

Now, these programs will return what are known as High Scoring Pairs (HSPs), the regions of your query sequence(s) that align well to the target. There are various options you can play with to improve the sensitivity or the specificity but a discussion of those would require far more details about what you're doing and would be best in a new question as well.

Once you have your HSPs, you can relatively easily parse them to select regions with a given range of sequence similarity values and of a specific length. Once again, that would be better discussed in a separate question once you have your results and can show an example.


How to get a list of genes corresponding to the list of SNPs (rs ids)? Is there a way for me to get a list of genes given a list of SNP rs ids? I found several questions asked with a similar goal years ago, and the answers are always about using multiple online tools, and/or R programming.

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

The relevant question I found (that are asked and answered more than five years ago.)

https://support.bioconductor.org/p/39278/
https://www.biostars.org/p/10437/
https://www.biostars.org/p/413/

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

A simpler solution, I don't know... but this is at least one Python solution using Biopython's ELink method via NCBI's Entrez E-utils.

The Biopython library is flexible enough (they have an in-depth tutorial worth reading) to modify the code below to to fit project-specific needs.  For example, the mapping between snp id and gene id is lost with this bulk request method.

from Bio import Entrez
Entrez.email = "youremail@gmail.com"

snp_ids = ["rs8179065",
           "rs12899845",
           "rs8179066",
           "rs80357332",
           "rs80356920"]

record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=",".join(snp_ids).replace('rs', ''),
                                  db="gene"))

This returns a record object

[{'LinkSetDb': [{'DbTo': 'gene', 'LinkName': 'snp_gene', 'Link': [{'Id': '374654'}, {'Id': '90381'}, {'Id': '672'}]}], 'ERROR': [], 'LinkSetDbHistory': [], 'DbFrom': 'snp', 'IdList': ['8179065', '8179066', '12899845', '80356920', '80357332']}]

which contains several iterable fields -- 'Link' has a list of the Gene UID's that can be converted to names using ESummary.

for gene_id in record[0]['LinkSetDb'][0]['Link']:
    handle = Entrez.esummary(db="gene", id=gene_id['Id'])
    uid_record = Entrez.read(handle)
    handle.close()

    uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
    print (gene_id['Id'], uid_summary['Name'])
    # ...or something more interesting than print

374654 KIF7
90381 TICRR
672 BRCA1

An option to preserve discrete mapping, is to loop over each snp id and extract the gene(s):

although I don't endorse the method below, especially for long lists of snps, you may flood the server with too many requests, maybe someone more familiar can comment on request limits or best practices for Entrez

for snp_id in snp_ids:
    record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=snp_id.replace('rs',''),
                                  db="gene"))
    results = record[0]['LinkSetDb'][0]['Link']
    for result in results:
        uid = result['Id']
        handle = Entrez.esummary(db="gene", id=uid)
        uid_record = Entrez.read(handle)
        handle.close()
        uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
        gene_name = uid_summary['Name']

        print(snp_id, uid, gene_name)

rs8179065 374654 KIF7
rs12899845 374654 KIF7
rs12899845 90381 TICRR
rs8179066 374654 KIF7
rs80357332 672 BRCA1
rs80356920 672 BRCA1


How is the GT field in a VCF file defined? As my question in SO was closed and asked to be posted in this forum, I am posting it here
I am not from the bioinformatics domain. However, for the sake of analysis, I am trying to pick up certain basics related to the GT field in the VCF file.
I know we have a field called Alleles. May I know under what circumstances GT takes a certain values and how they are called? Can you confirm my understanding?
Ref Alt    GT   Name
A    A     0/0  Homozygous
A    G     0/1  Heterozygous (does 1/0 also mean the same?) What's this 0 and 1 actually?
A   [C,CA] ??    ??
??   ??    1/1  HOM_ALT? How and why?

Can the experts here help me to fill the question marks and also help me understand with the Ref and Alt combinations when a genotype can take a value of 0/0 or 1/0 or 1/1 or 1/2 etc and what are the names for those values? Like when it is called home_alt etc
Any simple explanation for beginner like me (with no background in bioinformatics/biology) can be helpful

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

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


How to convert files to ADAM format? I would like to convert BAM and VCF files to ADAM format.

How do I do that?

These transformations are quite simple with the adam-submit script packaged with ADAM.

Transformation of BAM files

adam-submit -- transformAlignments sample.bam sample.alignments.adam

Transformation of VCF files

adam-submit -- transformVariants sample.vcf.gz sample.variants.adam


Using external list of PCs for clustering I am going to use principal components (PCs) comes from calcPCA function in URD program for clustering my cells in Seurat; So

seurat@dr$pca=URD_object@pca.load

> seurat <- FindClusters(
+     object = seurat,
+     reduction.type = "pca",
+     dims.use = 1:8,
+     resolution = 1.0,
+     print.output = 0,
+     save.SNN = TRUE
+ )
Error in GetDimReduction(object = object, reduction.type = reduction.type,  :
  cell.embeddings slot doesn't exist

There is seuratToURD function but there is not URD to seurat function I guess so I could use URD object for clustering in seurat
How I can overcome this error please? I want to see when I am using PCA from URD what would be my clustering by SNN algorithm in seurat

Thanks a lot for any help

There is a way to do this, and even better--there is documentation for how to do it! No surprise coming from the Satija Lab. In the vignette they perform multidimensional scaling, but the idea is the same. cmdscale() returns the cell embeddings. SetDimReduction() is the Seurat function you are looking for. No manual editing using @ required.

The authors use add mds embeddings, but if you want to add your own pca embeddings swap out "mds" for "pca" and mds for URD_object@pca.load.

# Before running MDS, we first calculate a distance matrix between all pairs
# of cells.  Here we use a simple euclidean distance metric on all genes,
# using object@scale.data as input
d <- dist(x = t(x = pbmc@scale.data))
# Run the MDS procedure, k determines the number of dimensions
mds <- cmdscale(d = d, k = 2)
# cmdscale returns the cell embeddings, we first label the columns to ensure
# downstream consistency
colnames(x = mds) <- paste0("MDS", 1:2)
# We will now store this as a new dimensional reduction called 'mds'
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "cell.embeddings",
    new.data = mds)
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "key",
    new.data = "MDS")

See the dimensional reduction vignette for more.


How to extract all pair end fastq.gz files from multiple subdirectory to a single folder I have all file names for each paired-end fastq pair as shown below. Could you please suggest how I can extract all fastq.gz files in a single folder?

sample name:

188833984]$ ls

Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc
Leaf_T1_FD_R2_L001-ds.7db8eb7e3426486db549426601b3a0bd
Leaf_T1_FD_R3_L001-ds.147177ecc03a46ccbbce3162d1185a0a

pair end file within sample name Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc

188833984]$ cd Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc

Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc]$ ls
Leaf-T1-FD-R10_S73_L001_R1_001.fastq.gz
Leaf-T1-FD-R10_S73_L001_R2_001.fastq.gz

I need all in fastq.gz format but the problem is that each sample, pair end R1 and R2 fastq reads in a single folder like this: and I need to extract all samples files in a single folder

Go to the common parent folder under which all FASTQ.GZ files reside. They could be in sub-folders under the folder, but all of them should have the folder you're in as a common parent folder. From such a folder, you should be able to run a find

find . -name "*.fastq.gz"

to get your FASTQ.GZ files. Now, there are surely better ways of doing this but here's one straightforward way.

Let's say you need to extract the FASTQs to ~/target/dir/, you'd do something like:

for f in $(find . -name "*.fastq.gz")
do
    tgt_file="~/target/dir/${f%.*}" #removed .gz from the name
    gzip -dc $f >$tgt_file #pigz -dc would be better if you have pigz
done

EDIT: OP has clarified that they wish to copy the files as such, not decompress them.

There are two ways to achieve this goal - you can either make physical copies of the files (thereby doubling the storage space you're currently using for them), or you can create symbolic links to the files (which are like shortcuts pointing to the actual files). I'll show you both ways here:

Let's assume the files are all under /fastq/files

find /fastq/files/ -name "*.fastq.gz" -exec cp {} ~/target/dir/ \; #This creates physical copies

cd ~/target/dir
find /fastq/files/ -name "*.fastq.gz" | xargs -I v_f ln -s v_f #These lines create symbolic links that point to the original files. Do not delete/move/rename the original files or these will become invalid.


Trying to divide number in fourth TSV column by a certain number according to its corresponding identifier in the first column I was wondering if anyone could help with my little issue in manipulating tsv files in order to calculate estimated coverage across many tsv files.

I have tsv files that look something like this:

ID        Info    Info2  x
M80428    info    info   1000
T83241    info    info   8400
E81292    info    info   10200

I also have another tsv file with all possible IDs and a corresponding number like so:

ID        y
M80428    27311
E11123    46531
P38163    14213
A34197    12892
T83241    37416
Q21345    27321
E81292    13101

x in the first tsv are estimated read counts and the y in the second is the genome length of each ID.

I would like to divide the "x" by the "y" in the second tsv that has the corresponding ID in order to get the estimated coverage of each hit. e.g. the first row would be 1000/27311, second row 8400/37416. I would then want this to be written to the fifth column of the tsv.

I am working in bash shell mostly here, and think that the way I would probably do this is with awk and by making the second tsv an associative array but I can't work out how exactly to do so. The python package dict seems like another option but my Python skills are practically zero.

Thanks for the help in advance!

Python dictionaries are indeed a good tool for that.

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

#!/usr/bin/env python3
import sys

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

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

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

sys.exit(0)

Testing it:

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


How to manipulate protein interaction network from String database in R? How can I manipulate protein-interaction network graph from the String database using STRINGdb bioconductor package and R?

I have downloaded the whole graph for Homo sapiens from STRING, which has about 20.000 proteins.

How do I read the file using that package?
How do I filter things I don't need? Supposing that I want to keep tumor data, as an example.

I think the easiest way is to download the graph using STRINGdb.

library(STRINGdb)
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )
full.graph <- string_db$get_graph()

Now you can use igraph, to manipulate the graph. Let's assume you want to take 200 proteins with the highest degree, i.e. number of edges they have.

library(igraph)

# see how many proteins do you have
vcount(full.graph)

# find top 200 proteins with the highest degree
top.degree.verticies <- names(tail(sort(degree(full.graph)), 200))

# extract the relevant subgraph
top.subgraph <- induced_subgraph(full.graph, top.degree.verticies)

# count the number of proteins in it
vcount(top.subgraph)

How to get disease specific genes?

There's no GO annotation for cancer or Alzheimer's disease. It is out of scope of the GO consortium.

What you can do, you can either take KEGG Pathways annotation, or manually select list of relevant GO-terms. Or acquire the list from one of the papers. For example annotation term 05200 corresponds to the cancer KEGG pathway. You can easily retrieve proteins associated with the annotation:

cancer.pathway.proteins <-
    string_db$get_term_proteins('05200')$STRING_id

And then perform subgraphing as described above.

Alternatively you can try to get an enrichment score for an every gene given it's neighbors (the way enrichment is shown on the string-db website). Then you can keep only those having top enrichment scores. Probably get_ppi_enrichment_full or get_ppi_enrichment functions will help you to do that.