LRT or LRT-like test on cyclical (Sleep) data I have RNA-Seq data from 4 time points (3 hours awake, 9 hours awake, 3 hours asleep, 9 hours asleep).  I'm interested in doing something similar to a LRT where genes are found to be significant if they changed at any point in the time course.

However, I'm not sure how to do this considering this is cyclical.  Running a LRT test in DESeq2 gives different results depending on which timepoint I choose as the beginning of the time course since I assume it doesn't calculate the last transition back to being awake (Or asleep depending on your reference).
Any ideas on an appropriate way to do this?

time = c(3,9,15,21)
ddsLRT <- DESeqDataSetFromTximport(txi.rsem.lrt, sampleTable, ~time)
ddsLRT<- DESeq(ddsLRT, test="LRT", reduced = ~1)

I think the problem is that you have time as a linear value rather than a factor. While this naively makes sense (after all, you have to go through 3 hours awake before you can get to 9 hours awake), what you end up doing is fitting a single coefficient to each gene of "linearly changes with time". What you want instead is factor, so you can have a system where things like per and clk can cycle:

time = factor(c(3,9,15,21))

The rest is then correct. Note that LRTs are not asking, "is there a change at any time point?", since there will likely be DE genes according to this that are not DE at any single time point. Rather an LRT such as this is asking, "What are genes that vary significantly over time?" There's a subtle different between those two questions. The former question is actually answered by a bunch of pairwise comparisons (taking the union of the DE genes).


How can I restrict the search space when querying GEO? I am using the following command to retrieve GEO Datasets entries for daf-2 c. elegans that have as their experiment type High Throughpout sequencing:
esearch -db gds -query "daf-2" | efilter -query "expression profiling by high throughput sequencing [DataSet Type] AND Caenorhabditis elegans [organism]"

However, among the results that I got, I noticed that I got this:
23. Impaired Insulin-/IGF1-Signaling Extends Life Span by Promoting Mitochondrial L-Proline Catabolism to Induce a Transient ROS-Signal
(Submitter supplied) Transcriptome profiling of three models with impaired insulin/IGF1 signaling. 1. Deep sequencing of endogenous mRNA from Caenorhabditis elegans N2 var. Bristol (wildtype) and daf-2(e1370) mutant; 2. Deep sequencing  of endogenous mRNA from murine embryonic fibroblasts (MEF)  wildtype and irs1-/- knockout; 3. Deep sequencing of endogenous mRNA from murine embryoinic fibroblast (MEF) insr+/- -lox and insr+/- knockout  Jena Centre for Systems Biology of Ageing - JenAge (www.jenage.de)
Organism:   Mus musculus; Caenorhabditis elegans
Type:       Expression profiling by high throughput sequencing
Platforms: GPL11002 GPL13776 14 Samples
FTP download: GEO (CSV) ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE36nnn/GSE36041/
SRA Run Selector: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA151765
Series      Accession: GSE36041 ID: 200036041

As you can see, in the Organism field, there are both Mus musculus and Caenorhabditis elegans. How can I restrict my search so that I only get entries that exclusively work on Caenorhabditis elegans? One naive way to do this would be to include NOT Mus Musculus in the query, but of course that would mean that I could get something else too.
Another way I can think of solving this is just to write a script that would do the extra filtering using regexes, but I was wondering if there's a simpler solution using the e-utilities functionalities

You can use the following query to first search GEO to get a list of GEO Series of interest, then find all linked SRA runs and perform a second filter step to keep only the SRA runs that satisfy another set of criteria as shown below.
esearch \
  -db gds \
  -query "daf-2[all fields] \
    AND expression profiling by high throughput sequencing [DataSet Type] \
    AND Caenorhabditis elegans [organism]" \
  | elink \
  -db gds -target sra \
  | esearch \
  -query "(#3) \
    AND Caenorhabditis elegans [organism]" \
  | efetch -format runinfo

Here, I am combining two independent queries using the (#3) term in the second query. More information about this is in the "Combining Independent Queries" subsection of "Searching and Filtering" section of the Entrez Direct documentation.


How can i modify FASTA headers in a multi fasta file using BioPython SeqIO I have a multi fasta file similar to this (relatively new here so uncertain of best way to present this; I have gone for an output and the code i used to make it - belt and braces...):

# create a test dataset
with open('test_dat.txt', 'w') as test_dat:
    l1 = [">AAA a aa", ">BBB b bb", ">CCC c cc"]
    l2 = ["ATGATG", "GACGAC", "TTGTTG"]

    for i in range(len(l1)):
        test_dat.write(l1[i] + "\n")
        test_dat.write(l2[i] + "\n")

# which outputs something like this
>AAA a aa
ATGATG
>BBB b bb
GACGAC
>CCC c cc
TTGTTG

I want to modify the fasta headers in my fasta file by adding some text to the beginning and, ideally, stripping out spaces (as i was advised that white space in the headings could cause problems for me down the line). My plan was to do this using BioPython and SeqIO.

I started like this:

# my attempt at altering the fasta headers using SeqIO
from Bio import SeqIO

# open my test data and create an output file
with open('test_dat.txt', 'r') as inputs, open('test_dat_out.txt', 'w') as outputs:

    # create dictionary of the fasta sequences using SeqIO
    my_dict = SeqIO.to_dict(SeqIO.parse(inputs, "fasta"))

    # renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header and remove white space
    for v in my_dict.values(): # using the values() to loop through...
        v.description = "TEST~~~TEST~~~" + v.description.replace(" ", "_")

    # write result to file
    SeqIO.write(my_dict.values(), outputs, 'fasta')

which gave me a result with some duplication in the fasta header:

#result i want
>TEST~~~TEST~~~AAA_a_aa
ATGATG

#result i get
>AAA TEST~~~TEST~~~AAA_a_aa
ATGATG

If i modify both the .id and .description in exactly the same way, as suggested here, but without removing white space:

# renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header
    for v in my_dict.values(): # using the values() to loop through...
        v.id = "TEST~~~TEST~~~" + v.id
        v.description = "TEST~~~TEST~~~" + v.description

I get almost what i want, but the header still has spaces:

>TEST~~~TEST~~~AAA a aa
ATGATG

I have tried to make the .id and .description match by adding a "_" to the end of .id, or by adding .replace(" ", "_") to both .id and .description (despite there being no spaces in .id), but I just cant seem to get exactly the output i want.

I would really appreciate any pointers on this.

You're very close, but unfortunately for legacy reasons this aspect of the Biopython parser is not very intuitive. The key point is the .id is anything up until the first white space character. One day I hope the .description and .id attributes are replaced by a single attribute like .header. Anyway, here is a solution for your case:

from Bio import SeqIO

to_add = "TEST~~~TEST~~~"

with open("out.fa", "w") as outputs:
    for r in SeqIO.parse("in.fa", "fasta"):
        r.id = (to_add + r.description).replace(" ", "_")
        r.description = r.id
        SeqIO.write(r, outputs, "fasta")

Output on your sample input:

>TEST~~~TEST~~~AAA_a_aa
ATGATG
>TEST~~~TEST~~~BBB_b_bb
GACGAC
>TEST~~~TEST~~~CCC_c_cc
TTGTTG


TCGA: gene IDs not appearing and other concerns I want to download RNA-seq datasets of 4 or 5 different types of cancer from the TCGA to investigate what is happening with my gene of interest. The problem is that I'm a physician and I don't know people used to deal with Big Data to help me.
I'm using a Bioconductor program called TCGAbiolinks and I'm able to download some datasets using the following script:
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("TCGAbiolinks")

browseVignettes("TCGAbiolinks")

library(TCGAbiolinks)
library(dplyr)
library(DT)

query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Gene expression",
                  file.type  = "normalized_results",
                  data.type = "Gene expression quantification",
                  experimental.strategy = "RNA-Seq",
                  barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
                  legacy = TRUE)

GDCdownload(query, method = "api", files.per.chunk = 10)

# GDCprepare transforms data in R format

# The default GDCprepare output is a SummarizedExperiment (SE) object,
# but you can also set to output the data frame.

data <- GDCprepare(query, summarizedExperiment = FALSE)

The problems is that: 1) the gene list is not appearing in the first column; 2) I don't know how many datasets I'll be able to download and analyse without a server.. I'm working with my Macbook..
Can someone familiar with the TCGA website or with the TCGAbiolinks help me with this project? I've been trying to progress in the last months, but the TCGA is a complex website and I'm feeling lost. I use R, I don't know SQL or other programming languages..

I had to play around a little to understand the object, but here's what I got:
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(TCGAbiolinks)
> query <- GDCquery(project = "TCGA-BRCA",
+                   data.category = "Gene expression",
+                   file.type  = "normalized_results",
+                   data.type = "Gene expression quantification",
+                   experimental.strategy = "RNA-Seq",
+                   barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
+                   legacy = TRUE)
--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg19
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-BRCA
--------------------
oo Filtering results
--------------------
ooo By experimental.strategy
ooo By data.type
ooo By file.type
ooo By barcode
----------------
oo Checking data
----------------
ooo Check if there are duplicated cases
ooo Check if there results for the query
-------------------
o Preparing output
-------------------

> GDCdownload(query, method = "api", files.per.chunk = 10)
Downloading data for project TCGA-BRCA
GDCdownload will download 3 files. A total of 1.306357 MB
Downloading chunk 1 of 1 (3 files, size = 1.306357 MB) as Wed_Aug_19_12_20_11_2020_0.tar.gz

Downloading: 580 kB

> data <- GDCprepare(query, summarizedExperiment = TRUE)
------------------
oo Reading 3 files
------------------
|====================================================|100%                      Completed after 0 s
------------------
oo Merging 3 files
------------------
Accessing grch37.ensembl.org to get gene information
Downloading genome information (try:0) Using: Human genes (GRCh37.p13)
Starting to add information to samples
 => Add clinical information to samples
 => Adding TCGA molecular information from marker papers
 => Information will have prefix 'paper_'
brca subtype information from:doi.org/10.1016/j.ccell.2018.03.014

> head(data)
class: RangedSummarizedExperiment
dim: 6 3
metadata(1): data_release
assays(1): normalized_count
rownames(6): A1BG A2M ... RP11-986E7.7 AADAC
rowData names(3): gene_id entrezgene ensembl_gene_id
colnames(3): TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07
  TCGA-AC-A23H-01A-11R-A157-07
colData names(83): barcode patient ... paper_PARADIGM Clusters
  paper_Pan-Gyn Clusters

I see from this representation that I can use metadata(data), assays(data), etc to get more information.
> head(SummarizedExperiment::assays(data))
List of length 1
names(1): normalized_count

Alright, so this is a named list. I know how to deal with lists: I can use either the index ([[1]]) or the name ($normalized_count)
> head(SummarizedExperiment::assays(data)$normalized_count)
             TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-01A-11R-A157-07
A1BG                             121.9879                      38.5795                      13.6874
A2M                             9703.8622                   37640.8373                    4767.0490
NAT1                            8488.0727                     114.8582                     270.4753
NAT2                               6.0583                       4.2938                       2.3804
RP11-986E7.7                   33016.6641                   11308.3460                    3894.9639
AADAC                              1.1359                     219.6976                       1.1902

I can now get more information on the samples using colData(). I'm truncating the output below because there are 83 columns.
> SummarizedExperiment::colData(data)
DataFrame with 3 rows and 83 columns
                                                  barcode      patient           sample shortLetterCode          definition sample_submitter_id         sample_type       state
                                              <character>  <character>      <character>     <character>         <character>         <character>         <character> <character>
TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046 TCGA-AN-A046-01A              TP Primary solid Tumor    TCGA-AN-A046-01A       Primary Tumor    released
TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-11A              NT Solid Tissue Normal    TCGA-AC-A23H-11A Solid Tissue Normal    released
TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-01A              TP Primary solid Tumor    TCGA-AC-A23H-01A       Primary Tumor    released
                                                        sample_id   is_ffpe sample_type_id  tissue_type submitter_id days_to_collection oct_embedded                pathology_report_uuid
                                                      <character> <logical>    <character>  <character>  <character>          <integer>  <character>                          <character>
TCGA-AN-A046-01A-21R-A034-07 57aff09f-0e97-4e65-a227-3dc7a8516367     FALSE             01 Not Reported TCGA-AN-A046                 27         true 1304FB17-A20A-4EBC-9CC8-1554808AC1F6
TCGA-AC-A23H-11A-12R-A157-07 7df59ca8-7e51-4581-9d7f-8bba0395ce17     FALSE             11 Not Reported TCGA-AC-A23H                478        false                                   NA
TCGA-AC-A23H-01A-11R-A157-07 d7e3b628-d5fd-4e79-9c4a-6409330fb8a7     FALSE             01 Not Reported TCGA-AC-A23H                478        false A7C7D409-D086-4A9B-8C8F-E7E231D5891D

I can get more information on the genes using rowData():
> SummarizedExperiment::rowData(data)
DataFrame with 19947 rows and 3 columns
                  gene_id entrezgene ensembl_gene_id
              <character>  <integer>     <character>
A1BG                 A1BG          1 ENSG00000121410
A2M                   A2M          2 ENSG00000175899
NAT1                 NAT1          9 ENSG00000171428
NAT2                 NAT2         10 ENSG00000156006
RP11-986E7.7 RP11-986E7.7         12 ENSG00000273259
...                   ...        ...             ...
RASAL2-AS1     RASAL2-AS1  100302401 ENSG00000224687
LINC00882       LINC00882  100302640 ENSG00000242759
FTX                   FTX  100302692 ENSG00000230590
TICAM2             TICAM2  100302736 ENSG00000243414
SLC25A5-AS1   SLC25A5-AS1  100303728 ENSG00000224281

You want the assay data as a data.matrix, I guess. That'd be:
> data_as_matrix <- data.matrix(SummarizedExperiment::assays(data)$normalized_count)

My session Info:
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] TCGAbiolinks_2.16.1

loaded via a namespace (and not attached):
  [1] bitops_1.0-6                matrixStats_0.56.0          bit64_0.9-7.1               doParallel_1.0.15           RColorBrewer_1.1-2          progress_1.2.2
  [7] httr_1.4.2                  GenomeInfoDb_1.24.2         backports_1.1.8             tools_4.0.2                 R6_2.4.1                    DBI_1.1.0
 [13] BiocGenerics_0.34.0         colorspace_1.4-1            tidyselect_1.1.0            gridExtra_2.3               prettyunits_1.1.1           bit_1.1-15.2
 [19] curl_4.3                    compiler_4.0.2              rvest_0.3.6                 Biobase_2.48.0              xml2_1.3.2                  DelayedArray_0.14.0
 [25] rtracklayer_1.48.0          scales_1.1.1                survMisc_0.5.5              readr_1.3.1                 genefilter_1.70.0           askpass_1.1
 [31] rappdirs_0.3.1              stringr_1.4.0               digest_0.6.25               Rsamtools_2.4.0             foreign_0.8-80              R.utils_2.9.2
 [37] rio_0.5.16                  XVector_0.28.0              pkgconfig_2.0.3             dbplyr_1.4.4                readxl_1.3.1                rlang_0.4.7
 [43] ggthemes_4.2.0              RSQLite_2.2.0               generics_0.0.2              zoo_1.8-8                   jsonlite_1.7.0              BiocParallel_1.22.0
 [49] zip_2.0.4                   car_3.0-8                   dplyr_1.0.0                 R.oo_1.23.0                 RCurl_1.98-1.2              magrittr_1.5
 [55] GenomeInfoDbData_1.2.3      Matrix_1.2-18               Rcpp_1.0.5                  munsell_0.5.0               S4Vectors_0.26.1            abind_1.4-5
 [61] lifecycle_0.2.0             R.methodsS3_1.8.0           stringi_1.4.6               carData_3.0-4               SummarizedExperiment_1.18.1 zlibbioc_1.34.0
 [67] plyr_1.8.6                  BiocFileCache_1.12.0        grid_4.0.2                  blob_1.2.1                  parallel_4.0.2              ggrepel_0.8.2
 [73] forcats_0.5.0               crayon_1.3.4                survminer_0.4.8             lattice_0.20-41             haven_2.3.1                 Biostrings_2.56.0
 [79] splines_4.0.2               GenomicFeatures_1.40.0      annotate_1.66.0             hms_0.5.3                   knitr_1.29                  pillar_1.4.6
 [85] ggpubr_0.4.0                GenomicRanges_1.40.0        ggsignif_0.6.0              codetools_0.2-16            biomaRt_2.44.1              stats4_4.0.2
 [91] XML_3.99-0.5                glue_1.4.1                  downloader_0.4              data.table_1.13.0           vctrs_0.3.2                 selectr_0.4-2
 [97] foreach_1.5.0               cellranger_1.1.0            gtable_0.3.0                openssl_1.4.2               purrr_0.3.4                 tidyr_1.1.0
[103] km.ci_0.5-2                 assertthat_0.2.1            ggplot2_3.3.2               openxlsx_4.1.5              xfun_0.16                   xtable_1.8-4
[109] broom_0.7.0                 rstatix_0.6.0               survival_3.1-12             tibble_3.0.3                iterators_1.0.12            KMsurv_0.1-5
[115] GenomicAlignments_1.24.0    AnnotationDbi_1.50.1        memoise_1.1.0               IRanges_2.22.2              ellipsis_0.3.1


Is the genome a programming language (i.e. LISP)? Can we analyze it with computer science theory? Question moved here from the biology stack exchange.

It looks like the genome is a kind of LISP language.

Operon structure of genome:
- https://github.com/philschatz/microbiology-book/blob/master/resources/OSC_Microbio_11_07_Operon.jpg

LISP function structure:
- http://support.ircam.fr/docs/om/om6-manual/res/listprefix.png

Nested genes:
- https://player.slideplayer.com/31/9782148/data/images/img11.jpg

Nested LISP cons cells:
- https://www.tutorialspoint.com/lisp/images/treestructure.jpg

The basic similarity I am seeing is of general structure: (func var0 var1 ... varn) and nested cons cells.

The correspondence with genome structure is the operon portion, which has:

operator = function
structural genes = variables
structural genes can be nested like cons cells

Is this a global structure for the genome, then a specific prediction is we can do something simple like counting open and close markers in the nested genes and they will be equal.

So, why would such a correspondence matter, if true?  The basic idea is that if the LISP structure is consistent, then we could maybe look for other programming constructs built on the structure, such as recursion/looping, maybe templating or even object oriented programming. Or just identify simple properties like balanced opening closing markers. The mere fact it has a LISP structure would be fascinating if true.  And since LISP is essentially an abstract syntax tree (AST), which all programming languages have, perhaps this structure I see in the genome is actually an AST, which implies an even more fundamental correspondence between the genome and programming languages.

More generally, we have a wealth of knowledge reverse engineering human programming languages.  If the genome is also a programming language similar in some respects to human programming languages, such as being a LISP variant, then we could apply this same reverse engineering knowledge base to better understand the genome.

So, first of all, is my high level observation correct?  Does the genome, at least for bacteria, have a global LISP like structure?  If yes, then if the genome is a kind of programming language, i.e. LISP, can we analyze it like software to understand its functionality?  Has anyone done this before?  Simple searching on Google Scholar turns up lots of software for analyzing genomes, but nothing about analyzing the genome as software.  It would be interesting to know whether anyone has seriously examined the question before, and either made some useful discoveries, or found out there is only a very loose analogy with programming languages, and thus knowledge regarding programming languages is not of use regarding the genome.

Hopefully this question is clear enough now, if not please ask me for further clarification.  I would like to get some sort of decent answer.

UPDATE: Someone who has put thought into "DNA is code" idea.  At a more general level than whether specific genome structures match specific programming structures.

UPDATE #2: The late Scott Federhen of NCBI wrote an article comparing genomic replication with the lambda operator.

UPDATE #3: It might be the case that LISP was inspired by the genome.

I see what the metaphor has been inspired by, where I think it got turned around, and might have a better one.

I'd like to back up a bit and establish some common ground.  It isn't wrong to expect that some aspects of genetic code and software code are similar, and while biology is usually very specific, some aspects can be seen productively from the light of computer science ideas.  For example, mass-action molecular dynamics, which are often a good model for protein interactions, are Turing-complete.  Here's a fun example where people have written a programming language out of molecular dynamics:  https://arxiv.org/abs/1809.07430

It's also true that proteins have particular substructures leading to certain functions, such as using alpha-helices to react less with the outside world.  Learning about larger substructures, called protein domains, help identify the functions of families of proteins.

The problem is that the parts of the genome which appear nested like LISP cells actually aren't included in the protein.  The regulatory regions, promoters, inhibitors, introns, and so forth never become part of the key "computational" machinery.  This nesting isn't actively part of the computational activity proteins do perform.  The main quality of LISP is that the structure and function are perfectly self-similar throughout an expression, but the structure of the regulatory regions are neither replicated down into the protein nor are necessarily recapitulated on a larger scale.

Unfortunately, we don't expect the computational structure of proteins to be the same from one protein to another.  This is because "binds only with X" and "binds only with Y" must be implemented differently to successfully be distinct operations, and it is these binding constraints by which molecular dynamics implement "computational" operations, such as "not X" and "not Y", and so those are distinct as well.

I think a better way to think about the regular structures of the regulatory apparatus is more like a network protocol.  Along an otherwise undifferentiated binary channel, one has to identify when and where the data is delivered to an application.  Similarly, along a strand of bases, the machinery of translation has to have an indication of where the "application" of the protein starts and ends, and when that protein should be delivered to the cell.  We might compare some of the images of genetic structure that provided in the original question with layouts of network headers: http://slideplayer.com/slide/4798836/15/images/16/TCP/IP+Packet+Structures.jpg  We might compare bacterial versus multi-cellular regulation as having different delivery mechanisms and data structures for delivering the application information, for different functional consequences, in the same way TCP and UDP have different structures for different trade-offs.

This is all complicated somewhat by the fact that proteins do then interact with the regulation of genes, leading to interesting "computational" dynamics such as feed-forward regulation.  However, I'm not aware of any particular structure of the regulatory apparatus that necessarily correlates with these larger functional patterns.  The opposite is frequently biologically useful: the same structure of a protein helping to activate one gene might simultaneously inhibit another.

One book I've found very accessible and useful in understanding where computation is useful biologically is Dennis Bray's "Wetware: a computer in every living cell".  It gave me a clearer picture of how computation occurs in cells.  That should open up some introductions to systems biology which explains the biological-role of larger molecular "computational" patterns.


How to use analysis output in my pipeline I am running an analysis and I want to edit the output that I received from that analysis. However I cannot define the output name in my pipeline, because after running the analysis it is asking me to name the output file. How can I include that output file in my pipeline ?

I assume you are referring to the pepstats program from EMBOSS. The program has a -outfile flag that allows you to specify your desired name for the output file when you're not running the program interactively.


Extracting specific sections from FASTQ determined by positions Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions.

For example, say I want to extract the nucleotides at position 1, 3, 5-7 from each sequence:

Input:

@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789

Output:

@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

Note, the phred-scores are constructed such that they correspond to the position for easier reading.
Note2, the phred-score subsections should be included accordingly (i.e. the same positions)

Ideally, a command should work this way:

subselect 1,3,5-7 input.fq > output.fq

Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5` and pasting) but maybe out of an box solution does already exist which evaded me.

For the specific regions you mentioned, you could do this in awk. It won't be fast for a large file, but it will work:

$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

$0 is the current line. The script above will replace the line with three substrings of itself, one starting at position 1 with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The NR%4==2 || NR%4==0 ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.

IMPORTANT: this assumes you always have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine.

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (1,3,5-7):

#!/usr/bin/env perl

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

Save that as foo.pl and then run like this:

perl foo.pl file.fq 1,3,5-7


How to convert a .fastq chip-seq file into a wig/bigwig track and how to extract chip-seq peaks? I have a reference .fasta file and a raw .fastq file with chip-seq data. I am trying to create a bigwig track from the and .fastq and .fasta ref file of the raw signal. Then I would like to do some peak analysis for the chip-seq track.

Can someone help me with googling the write things, I feel really lost while trying to google it and not sure where to start from. I was thinking to do .fastq -> .bam -> wig -> bigwig and then somehow to do peak detection analysis.

Before you proceed, it is usual to do quality control on the .fastq data.  The FASTQC program is a good place to start.

Once you have checked the quality, trimming is the typical next step to remove low-quality data.  There are options such as Trimomatic, trimGalore for this.

Once you are sure you have clean data, your first step is alignment - align the sequences in the .fastq file to the reference sequence in the .fasta file.  bowtie2 is one option for this process.  This gets you to the .bam stage.

These steps will take a little time, some research, and may be somewhat frustrating, but they are the foundation for what you want to do. One option for somewhat automating these steps is to use the Galaxy platform (usegalaxy.org), otherwise these are mostly command-line tools.

The bam to wig and wig to bigwig steps are relatively simple conversions, and Biostars has answers for these.


Prediction of prokaryotic origins of replication (ORI) I want to predict origins of replication (ORI) on hundreds of prokaryotic genomes. The most straight-forward solution would be to use most commonly used tool, Ori-Finder.

It uses integrated gene prediction, analysis of base composition asymmetry, distribution of DnaA boxes, occurrence of genes frequently close to oriC regions and phylogenetic relationships to predict ORIs.

However, I came across a paper by Parikh et al which describes a single measure (correlated entropy measure / CEM) with which they achieved 100% ORI-identification on 500 genomes.

I like the CEM approach more than a close-source Windows-only black-box.

Parikh tested their algorithm on NCBI genomes. How are ORI annotated on these reference genomes? I couldn't find such an annotation. Does the sequence simply begin at the ORI in leading strand direction by convention? Were the ORI of the reference genomes identified by experiment?

I have not seen the papers but I know the technique SVM. This is a machine learning approach its pretty good, but has been around this the 1980s.

Its quite simple known ORIs are given with the associated genetic trait, e.g. GC skew, SVM will then regress the GC skew of a known ORI against the GC content of the genome, there will be a training, parameter optimization and test split in the data, typically 80:20:20. De novo data will then be feed into the trained algorithm, in the same way the test data was, and a prediction will be generated. There are weaknesses, the authors have used only one ML technique, typically numerous will be used e.g. SVM, random forest, lasso regression etc... In addition, the prediction accuracy for GC is low (75%), in addition 500 genomes isn't great for ML, you need 1000s.


how to obtain GRCh38 bed and symbol files for HOCOMOCOv11 I am trying to a recently published software called diffTF and contacted the authors on how to obtain the input HOCOMOCO files for GRCh38. While I wait for the answer, I thought I would post the question here as well:

The tool requires bed files for each of the transcription factors (TF) in the genome, such as this bed file (TAB separated), which seems to give the locations of different TF binding sites in the genome:

chr1    3457132     3457141     GTTGCGTGA   1266   +
chr1    4595374     4595383     GTTGCGTGA   1266   -
chr1    4996631     4996640     TTTGCGTGC   1276   +
chr1    5005357     5005366     TTTGCGTGC   1276   +
chr1    5170456     5170465     TTTGCGTGC   1276   -
chr1    5567302     5567311     GTTGCGTGA   1266   -
chr1    8676279     8676288     TTTGCGTGC   1276   -
chr1    8859719     8859728     GTTGCGTGA   1266   +
chr1    9600920     9600929     GTTGCGTGA   1266   +
chr1    9650971     9650980     GTTGCGTGA   1266   +
...

Then there is a symbols file to map the TF to their gene(?) that looks like this:

SYMBOL  ENSEMBL HOCOID
Ahr     ENSMUSG00000019256      AHR
Aire    ENSMUSG00000000731      AIRE
Alx1    ENSMUSG00000036602      ALX1
Ar      ENSMUSG00000046532      ANDR
Tfap2a  ENSMUSG00000021359      AP2A
Tfap2b  ENSMUSG00000025927      AP2B
Tfap2c  ENSMUSG00000028640      AP2C
Tfap2d  ENSMUSG00000042596      AP2D

I googled for HOCOMOCO and found a bunch of pwn or words files, but none of them seems to be similar to the bed file above. Any ideas where to get these for GRCh38?

I guess you're right that you need predictions of transcription factor binding site (TFBS) through genome.
There're two options to get them:

You can generate predictions using any motif scanning tool. E.g. you can use our tool SPRY-SARUS which is able to generate bed-file in necessary format.
In this case, you can choose more or less strict thresholds for each prediction (i.e. more false positives vs more false negatives). Standard thresholds for HOCOMOCO motifs can be found at its downloads page (we, as HOCOMOCO authors, recommend thresholds corresponding to P-value of 0.0005).
I had no experience with diffTF but it seems that authors prepared TFBS predictions and put them on their site. Please, take a look at this page of manual. If you don't need extra-flexibility in the generation of TFBS predictions, it can be the simplest way to proceed.


Off-by-one error when mapping Ensembl IDs to gene symbols I have a matrix with ensemblIds which I am trying to convert to gene symbols using BioMart. Some of the ensemblIds will not be mapped, so I need to subset the rows of the matrix. It seems like BioMart adds one ensemblId or at least I see one more result ultimately than should be. Here is what I am trying to do:

library(dplyr)
library(tibble)
ids <- rownames(mtx_combined)
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id",
                           "mgi_symbol","description"),values=ids,mart= mart)
mtx_mapped <- as.matrix(as.data.frame(mtx_combined) %>%
rownames_to_column("gene") %>%
filter(rownames(mtx_combined) %in% list$ensembl_gene_id) %>%
column_to_rownames("gene"))

rownames(mtx_mapped) <- list$mgi_symbol   //Fails

Length of the returned list is 18337, number of rows in the initial matrix is 18376. After subsetting I expect to get 18337 rows, but getting one less - 18336. What am I doing wrong here?

There’s no 1:1 mapping between the IDs of different database schemas. As a consequence, some Ensembl gene IDs map to multiple MGI symbols, or to none (and vice versa).

Therefore, you can’t assume that the number of result rows will be equal to the number of query IDs.


Identification of Microorganism I have the 16S DNA Sequence of a Microorganism. Can I confirm its identity by only running a BLAST or do I have to follow other methods. If so please mention the methods.

As always with bioinformatics, any reasonable answer will be dependendant on what you want to do with this sequence, and you haven't provided enough information in that regard.

For most exploratory analysis BLAST is good enough. You don't have to follow any existing methods. The benefit of using NCBI's web BLAST server is that it provides you with a very quick way to find a match from the sequence information in GenBank.

BLAST does a local alignment. If you have the full 16s sequence, and it's a novel organism, you might get different match results from a program that does global alignment, but then again you might get different results from a different local alignment tool.


GRCm37-designed exon target enrichment, which reference to use? I have exomes from 24 individual mice. The exomes are the product of Roche's SeqCap EZ HyperCap target enrichment kit. I see that the mouse exome design comes from mm9/NCBI37, the previous major reference genome released in July 2007.

Which reference should I use to align these reads? What artifacts may arise from mapping these reads designed from mm9 to the mm10 reference? I have already aligned the reads to mm10 and noticed an erroneous exon right in the middle of an intron, perhaps a target that was an exon in the mm9 release but revised out of the mm10 release. My goal is to investigate structural variation.

The only possible issue with using GRCm38 as the reference is that since GRCm37 was used in the bait design there are likely a few probes that will show artificial structural variations, due to errors in GRCm37. However, these will be present in all of your mice and if you use liftOver to look at the regions around the probes then such regions should become obvious. In fact, that'd be a convenient thing to do if you want to keep your (likely better) GRCm38 alignments, flank each probe site by a few hundred bases and then use liftOver to get the GRCm37 coordinates. If you go through the result you'll be able to find the false-positive structural variants due to the baits being designed with an older reference.

Alternatively, just use GRCm37 and hope that any issues with it don't result in false-negative SV calls.


Preparing binary matrix data for Scikit classification algorithms I made this post in regular stack overflow but I was told about this awesome feature by @nbryans.

I am a researcher (my programming knowledge is small) conducting analysis on a set of antibiotic (methicillin) resistant and a set of antibiotic (methicillin) susceptible. The sole difference between these two sets is assumed to be the resistance (R) vs susceptibility (S). I want to find if a genetic element in the resistant or susceptible genomes is correlated to resistance (positively or negatively).

I have done a complete fisher exact test already and will detail how I did so below. However I would like to figure out how to prepare my data for a more complex learning algorithm such as those offered by sci-kit.

(many of my techniques were learned from this article)

First I used KMC to split each fasta file for both R and S genomes into k-mers (I don't have enough "reputation" so I can't use >2 links but k-mers are essentially genetic sequences of k length). The output for KMC was a dump of k-mers present in a genome.

Example: kmc dump

AAAAAAAAAACGAATGTACACAATCGAA    1
AAAAAAAAACGAATGTACACAATCGAAC    1
AAAAAAAAATCAAATCCTGACTATTTAG    1
AAAAAAAAATCGGTCAATTCATTAAAAG    1
AAAAAAAACAAACATGAAGACCTTGTTA    1
AAAAAAAACGAATGTACACAATCGAACA    1
AAAAAAAACGTGTTAAAGTGAATCACAC    1

I then merged the data into one binary matrix which had columns as genomes and rows as k-mers.

Example: 1 resistant genome and 1 susceptible genome

kmer    R_1    S_1
AAAAAAAAAACGAATGTACACAATCGAA    1   0
AAAAAAAAACGAATGTACACAATCGAAC    1   0
AAAAAAAAATCAAATCCTGACTATTTAG    1   0
AAAAAAAAATCGGTCAATTCATTAAAAG    1   0
AAAAAAAAATTCCCTTCTAATCTTGAAT    0   1
AAAAAAAACAAAAATTATATAAAGCGAA    0   1
AAAAAAAACAAACATGAAGACCTTGTTA    1   1
AAAAAAAACAACCACCCATACATTGAGT    0   1
AAAAAAAACCCTTACAACAAATATGTAA    0   1

This is the data I have been working with and would like to analyze further with sci kit or other classification algorithms. Essentially each k-mer is a feature whose correlation with resistance is based on its presence in R vs S genomes. This is where I would like help transitioning this data into a dataset format suitable for sci kit.

I conducted a Fisher Exact test on each k-mer using R. For each k-mer a 2x2 matrix was created with the first column as resistant and the second column as susceptible. The two rows were present or not. So four numbers (# times k-mer was present in R, # in susceptible, # not in R, # not in S).

      R    S
yes   #    #
no    #    #

I used the following code in R:

phenotype = as.numeric(grepl('^R_', colnames(raw_kmer_table)[2:ncol(raw_kmer_table)]))
fe_results = apply(kmer_pres_abs_matrix, 1,
               FUN = function(row) {
                 fe_mat = matrix(0, ncol = 2, nrow = 2)
                 fe_mat[1,1] = sum(row == 1 & phenotype == 1)
                 fe_mat[1,2] = sum(row == 1 & phenotype == 0)
                 fe_mat[2,1] = sum(row == 0 & phenotype == 1)
                 fe_mat[2,2] = sum(row == 0 & phenotype == 0)

                 fe = fisher.test(fe_mat)
                 return(fe$p.value)
               }
)

Now that you have read how my data is currently formatted I would like to hear suggestions on how I can fit this data into a more complex test using sci kit or another resource. Algorithms I am interested in using are any you see fit and previously field-tested algorithms such as adaboost and forest learning algorithms. Also please ask me any clarifications or things I may have left out of this post!

Since you are using R, you probably don't want to use scikit-learn, which is for Python. However, there is a similar R library mlr ("R package to make machine learning in R easy") that provides a unified interface to all popular machine learning methods. Check their tutorial on how to get started.


chromPlot R package: remove scales I am using chromPlot R package to plot some SNPs and their associated score, however, I cannot find a way to remove the scale of such score (circled in the attached picture). Does anybody know how to do it?

Thank you so much.

Setting the legChrom argument to NA will omit plotting a legend.

Found in the manual for chromPlot on Bioconductor for release 3.7.


Concatenate multiple fasta files into one file per unique entry I have eight fasta files, one per gene, with 106 fasta entries per file, like this:

gene1.fasta:

>sample1_434564
GATCTCGTAGATGTGA
>sample2_2346764
GATCTCGTAGATGTGA

gene2.fasta:

>sample1_532345
ATGCTGTGAAGTA
>sample2_6040103
ATGCTGTGAAGTA

What I want to do is create one fasta file per sample so that I get all gene sequences for that sample in one file. As you can see, the sample ID is the same before the _ in the fasta headers, but the rest is unique for each fasta file.

Is  there any way to do this with awk or other command-line tools?

Desired output:

Sample1.fasta

>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

Sample2.fasta

>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

You can do this simply enough in awk:

$ awk -F'_' '{
               split(FILENAME,n,/\./);
               if(/^>/){
                  sample=$1;
                  sub(/>/,"",sample);
                  $0=">"n[1]
               }
               print > sample".fasta"
             }' gene1.fasta gene2.fasta ... geneN.fasta

If I run this on the two files you gave, I get:

$ cat sample1.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

$ cat sample2.fasta
>gene1
GATCTCGTAGATGTGA
>gene2
ATGCTGTGAAGTA

Explanation

-F'_' : set the field separator to _ as an easy hack that lets us use $1 for the sample name (sample_NNN, so sample will be $1 and NNN will be $2).
split(FILENAME,n,/\./); : FILENAME always holds the name of the file currently being processed. Here, we are splitting it on . to get the gene's name. This assumes all of your files are named $geneName.fasta. The split creates the array n whose 1st element, n[1] is the gene name.
if(/^>/){ : if this is a fasta header line.
sample=$1 : set the variable sample to the sample's name. See the first note about -F'_'.
sub(/>/,"",sample); : remove the > from the sample variable. We don't want that in the file name.
$0=">"n[1] : set the current line to be > and the gene name.
print > sample".fasta" : print each (modified) line into the appropriate file name.


How to convert a list of genes to one element per row? I got two lists of about 300 genes from my GO analysis, like the following:

Mtss1/Sox11/Fbn2/Mmp2/Tcf7l2/Six2/Ahr/Fat4/Lef1/Stc1/Acvr2a/Glis1/Nog/Igsf10/Twist2/Igf1/...

How can I convert the list to one element per row (for Venn diagram  on Venny 2.1), as:

Mtss1
Sox11
Fbn2
Mmp2
...

$ echo "Mtss1/Sox11/Fbn2/Mmp2/..." | tr '/' '\n' > list.txt

For more information on tr: https://en.wikipedia.org/wiki/Tr_(Unix)


How to transfer gff annotations in genome with extensive duplications? Microbial genomes can contain extensive duplications. Often we'd like to transfer annotations from an annotated species to one that is newly sequenced.

Existing tools (e.g. RATT, LiftOver, Kraken) either make specific assumptions about how closely related the species are or fail to transfer when multiple matches are found in the new genome, especially if the sequences are highly similar.

Specifically, I have a synthetic biology application where genes can duplicate extensively. They are identical in sequence but duplicated many times and be relocated (i.e., not just adjacent to each other). None of the above mentioned tools are able to transfer coordinates of annotations to genomes with multiple copies of features.

Are there any pre-existing tools or software that transfer annotations in this scenario? Ideas for ways to do this robustly?

There is one very simplistic way I use which might work for what you are doing, it is similar to what terdon proposed.

Take a de-novo microbial genome annotation tool (I have my own, but you could use/modify prokka). Tools like these often first predict gene boundaries (with other tools like prodigal or glimmer) and then try to assign a function to found genes. This function assignment is often done with BLAST and other tools ... and that is where you can go in and modify to do what you need.

I use a "knowledge" protein database of genes I want to have very strictly annotated as a first line of annotation (e.g. in your case: the annotated genomes). For that I loop through very strict identity/similarity parameters which get gradually relaxed.

E.g.:
Loop 0: only transfer annotations at 100% DNA identity, same length.
Loop 1: only transfer annotations at 100% similarity, same length.
Loop 2: only transfer annotations at 99% similarity, length +/- 1%.
...
Loop n: only transfer annotations at 100-(n-1)% similarity, length +/- (n-1)%.

In each loop, obviously only annotate what has not been annotated in previous loops.

After that, use "normal" annotation pipeline of the tool to annotate the rest.


How to filter for "CAF=0.5747,0.4253" In my gzipped VCF file I have many variants, and they report allele frequency such as:

##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

and

CAF=0.5747,0.4253

I would like to filter out the variants >= 1% allelle frequency. How do I do that?

using vcffilterjdk:

you want to keep the variant if there is no INFO/CAF.

Otherwise, get All values of CAF as a String list (because it is declared as a string in the VCF header...), stream those data, skip the first value for REF , convert each string to double, you want to keep the variant if any allele is rare (< 0.01)

java -jar dist/vcffilterjdk.jar -e 'return !variant.hasAttribute("CAF") || variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf


How do I download the mitochondrial haplogroup datasets for human genetics online? I seem to have landed on the mitomap.org site, but I don't know what to make of it or what do with it / how to get the genomes onto my computer. It sounds like the genomes are stored in GenBank, but that mitomap simply lists which of the GenBank genomes are human mitochondrial haplogroup genomes.
First of all, I can't find a simple list of links or list of names of the haplogroups the mitomap has... Can you show me where I can get or download that? They say:

This brings our total number of FL sequences to 51,673, and the number of CR sequences to 74,660. Our SNVs now total 19,227.

I don't know what FL or CR (Control Region sequences?) means, but I am expecting to see lists of genomes or something. What am I to do with this information? Where do I find this data? Is this it? If I click on one of the results, I see this. It has path nuccore/MT742594.1, can I find this in some NCBI FTP server somewhere? I don't see anything of that folder structure in the GenBank FTP server...
Basically my question is:

Where is the complete list of human mitochondrial haplogroup genomes?
Where can I download all of those human mitochondrial haplogroup genomes so I can get them onto my computer?

Hmmm, I am not experienced enough in biology or bioinformatics to truly understand the question, but I am guessing you want to download the FASTA file or genebank file containing the sequence of the accession number you listed above.
There are 2 ways to do this -

Using your browser, if you want to download a small sample of files. On the right hand side of your screen you can see the "send to" button on this link, click on it and choose the what's appropriate for you and then click on create file.

https://www.ncbi.nlm.nih.gov/nuccore/MT742594.1

If you familiar with Biopython and python programming you can download multiple files at once in whatever format you want, there is an easy script to do this and has been taken from Biopython's documentation on biopython.org
import os
from Bio import SeqIO
from Bio import Entrez
Entrez.email = "A.N.Other@example.com"  # Always tell NCBI who you are
filename = "MT742594.1.gb"
if not os.path.isfile(filename):
 # Downloading...

 net_handle = Entrez.efetch(

     db="nucleotide", id="MG762674", rettype="gb", retmode="text"

 )
 out_handle = open(filename, "w")

 out_handle.write(net_handle.read())

 out_handle.close()

 net_handle.close()

 print("Saved")

print("Parsing...")
record = SeqIO.read(filename, "gb")
print(record)

You can make a list of the accession numbers of the sequences you want to download and pass the list in the id, if you are curious as to how the code works visit biopython.org .


Why are my genes filtered for Gene Ontology term enrichment? When I enter a list of gene names into common gene ontology (GO) term enricher (e.g. MsigDB or ToppGene) only a small fraction of genes is actually used for the enrichment:

I am struggling to understand why precisely only such a small number of genes is mapped here.

My dataset is mapped with the GENCODE annotation. What annotation is the base of gene ontology? (The picture mentions entrez genes)

What is "click here for details" saying? That link should contain an explanation of what happened.

GSEA usually requires entrez_id to run. If your list contains gene_symbols, then before enrichment, you need to convert the genes names. This can and is (always) tricky as genes are known with different names. For example, by googling I couldn't find any AASS3 gene, but I found AASS, which is also named LKRSDH,LORSDH, and LKR/SDH. Depending on the software you are using for conversion, and the genome version used some genes will be missed. I usually run (in R) a combination of geneSynonym and biomaRt to fill the gaps in the lists.

Even better, would be to use ensembl_id or even better entrez_id directly


Add tags to the read pairs in a bam file I have a bam file from paired-end sequencing, in which each R1 has the RX tag assigned (via UMItools group) denoting a Unique Molecular Identifier (UMI).
I want to call molecular consensus reads using fgbio but fgbio GroupReadsByUmi gives me an error, because it expects both R1 and R2 to have the RX tag set.
I would like a tool that will go through my bam file and assign each R2 the same RX tag as the R1.

The best thing would be to figure out how to do this in UMItools. I don't know that tool at all, but it seems that it is supposed to be able to deal with / tag both reads in paired data. e.g. did you use --paired flag, etc.
I would suggest investigating that route first before you look at the next part of my answer.

Per Devon Ryan's comment suggestion I would suggest a workflow similar to this:

sort bam file by query name (samtools sort -n) to get the R1s immediately following the R2s
Write a short script/command that

reads through the bam file and for each line
looks up the tag in the R1, prints R1 line
reads the R2, substitutes the RX tag in, prints R2 line

pipe that through samtools view and possibly samtools sort back to a bam file of the same type that you started with.

My pysam is a little rusty but something like this might work
import pysam

bam = pysam.AlignmentFile("bam", "rb")
# print header
print(bam.header)

tag = None
for read in bam:
    r2 = False
    if tag:
        read.assign_umi_tag() # assign to read- not sure where/how
        r2 = True
    else:
        tag = read.lookup_umi_tag()  # i don't know where this tag is in bam format
    print(read)
    if r2:
        tag = None

I haven't tested that in any way, caveat emptor.
I am not personally familiar with UMItools and how it assigns tags so can't suggest anything more specific, sorry.


Formula for k-mer coverage Let $C$ be base coverage, $R$ is the length of reads and $K$ is the length of $k$-mer. Then $k$-mer coverage $C_k$ can be computed as $C_k = C\cdot(R - K + 1)/R$.

Could someone please explain why is this equation valid (I'm mostly confused as why it is divided by $R$)?

Source: Velvet manual Section 5.1

I was still puzzled from the answers, so I tried to calculate with all the steps. I take this definition "$C_k$ is the number of reads containing a k-mer." and corresponding definition for coverage ($C$): "$C$ is the number of reads covering a base".

Coverage is $C = \frac{T \cdot R}{L}$, where $T$ is total number of reads, $R$ is read length and $L$ is length of genome. Given the $C_k$ definition, $C_k = \frac{T (R - K + 1)}{L-K+1}$, where $R - K + 1$ is just number of kmers in a read, and $L-K+1$ is number of kmers in a genome. Then,

$$C_k = \frac{T (R - K + 1)}{L-K+1} = \frac{T (R - K + 1)}{L-K+1} \cdot \frac{R}{R} = \frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L - K + 1}$$

since $L >> K$, we can approximate $L - K + 1 \approx L$, then we reduce the expression to

$$\frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L} = \frac{R - K + 1}{R} \cdot C$$

which is the formula for $C_k$.


What is the meaning of these misaligned reads in a sequencing run? I am analyzing some SARS-CoV-2 sequencing runs abd often find read alignments like the one in the image.
GACCCGGTGGGTTTTAC---ACTTAAAAACACAGTCTGTAC
GACCCTGTGGGTTTTAC---ACTTAAAAACACAGTCTGTAC
          GTTTCAC---ACTTAAAAACACAGTCTGTAC
          GGTTTAC---ACTTAAAAGCACAGTCTGTAC
CCACACTCTCCTAGCAC---ACTTAAAAACACAGTCTGTAC
CCACACTCTCCTAGCACCATACTTAAAAACACAGTCTGTAC

It basically says there are reads that share part of a sequence and part of another one. The first 4 reads are the same sequence and it's a sequence that you can effectively find in the official SARS-CoV-2 assembly but the other two reads are a mix of this one with another sequence (and both are the same mix) but with only a small (3bps) indel of difference.
The sequence of the last two reads does not appear anywhere in the official SARS-CoV-2 sequence.

What does it mean?
Are these sequencing errors?
If that is the case how is it that the same mistake is repeated with similar but not exactly the same data?
If these are real RNA sequences present in the probe, why are these not used anywhere in the final assembly and how did they arrive there?

Is there any other option that I am missing?

The unaligned parts of the bottom two sequences are from elsewhere in the sars-cov-2 genome if you blast them. Likely causes of this are library prep. artifacts (not terribly common, but you'll find a few cases here in there once you have millions of reads) where fragments end up ligating together or (probably more likely) rearrangements of the viral sequence. My guess would be that the latter is what occurred, given how rapidly viruses reproduce. Whether that arrangement can lead to viable viruses I don't know (that's too far out of my field).

Assemblies represent at best a snapshot of a single example of a species at a single moment in time. Actual individuals in the wild can vary quite substantially, especially when they don't need to be viable.


about the cpm normalization after using normalization factor Is it okay to use CPM normalization (with/without log transform) after using TMM normalization? Why do we need both?
 library(edgeR)
 library(SummarizedExperiment)
 load(url("http://duffel.rail.bio/recount/SRP049355/rse_gene.Rdata"))
  counts <- assays(rse_gene)$counts
  y <- as.matrix((counts))
 y <- DGEList(counts = y, group=c(1,2,3,4,5,6,7,8,9,10))
  y <- calcNormFactors(y)
 z <- cpm(y, normalized.lib.size=TRUE)
scaledata <- t(scale(t(z))) # Centers and scales data.
hc <- hclust(as.dist(1-cor(scaledata, method="spearman")), method="complete") # Clusters columns by Spearman correlation.
 TreeC = as.dendrogram(hc, method="average")
 plot(TreeC,
 main = "Sample Clustering",
 ylab = "Height")

In fact the cpm function will use the size factors (TMM) that were calculated with calcNormFactors. In the absence of the size factors cpm would do a naive per-million scaling so only correct for differences in library size. The power of the size factors is to also correct for library composition. I recommend this video for details. In short: What you do is fine. I would suggest though that you use the log = TRUE option for cpm as one typically performs downstream applications (even when using the Z-score) on the log scale. Or alternatively t(scale(t(log2(z+1)))).


VEP output SIFT_score unclear We have been experimenting with VEP (Variant Effect Predictor). One of the meta data attributes that we are interested in is the SIFT score, indeed when we apply the dbNSFP plug we get a column containing the scores (named SIFT_score). However, I don't understand why there are sometimes dots or multiple values in the fields. For example, the gene ENSG00000196924 below has 5 transcripts:

The SIFT_score column contains several values, not 1 per transcript/rs-number...

Here is another example that confuses me (I added the SIFT_pred column this time):

There are two mutations, the lower one can be expressed in 4 transcripts and I thus understand that there can be 4 SIFT scores, but why are all for in given in every row? Is the first one the SIFT_score for the first transcript?

One last example, again 4 transcripts, but now 2 of the scores are dots, what does that mean?

I have been looking for quite some time now how to interpret this data, any help is appreciated.

The dbNSFP plugin from VEP accesses tables of data for each variant from dbNSFP and pulls out the values. dbNSFP provide their SIFT scores in that format: a score for every transcript affected by the variant, all on one line. The lookup is just for the variant, not for the variant/transcript combo, so they provide scores for all variant/transcript combos. You can also get a column that gives you a list of the transcripts or proteins (Ensembl_transcriptid or Ensembl_proteinid) in order so you know which score goes with which transcript.

A better way to get SIFT scores with VEP is to get them directly from VEP, rather than using dbNSFP. This will get the SIFT score that goes with the transcript on the line with the relevant transcript.


What methods are available to find a cutoff value for non-expressed genes in RNA-seq? I have a gene expression count matrix produced from bulk RNA-seq data. I'd like to find genes that were not expressed in a group of samples and were expressed in another group.

The problem of course is that not all effectively non-expressed genes will have 0 counts due to sequencing errors, or because they were expressed in a small subset of cells.

I'm interested in solutions using R.

I'd like to find genes that were not expressed in a group of samples and were expressed in another group.

This is, fundamentally, a differential expression analysis, with a twist. To solve this, you’d first use a differential expression library of your choice (e.g. DESeq2) and perform a one-tailed test of differential expression.

Briefly, you’d perform the normal setup and then use

results(dds, altHypothesis = 'greater')

To perform a one-tailed test. This will give you only those genes that are significantly upregulated in one group. Check chapter 3.9 of the vignette for details.

Of course this won’t tell you that the genes are unexpressed in the other group. Unfortunately I don’t know of a good value to threshold the results; I would start by plotting a histogram of the (variance stabilised) expression values in your first group, and then visually choose an expression threshold that cleanly separates genes that are clearly expressed from zeros:

vst_counts = assay(vst(dds))
dens = density(vst_counts[, replicate])
plot(dens, log = 'y')

(This merges the replicates in the group, which should be fine.)

Counts follow a multimodal distribution, with one mode for unexpressed and one or more for expressed genes. The expression threshold can be set somewhere between the clearly unexpressed and expressed peaks:

Here I used identify(dens) to identify the threshold interactively but you could also use an analytical method:

threshold = identify(dens)
quantile = sum(dens$x < dens$x[threshold]) / length(dens$x)

# Using just one replicate here; more robust would be to use a mean value.
nonzero_counts = counts(dds, normalized = TRUE)[, replicates[1]]
nonzero_counts = nonzero_counts[nonzero_counts > 0]

(expression_threshold = quantile(nonzero_counts, probs = quantile))

26.5625%
4.112033


How to obtain file offsets corresponding to range query in block-gzipped file How can we retrieve the file offsets of the zipped-chunks in a block-gzipped file that contain records for a region/window of interest (e.g. Chr01:100-Chr02:145)?

It seems like the CLI of samtools, or tabix should support that, but I can't find a way to just get the offsets in the datafile, without having the accompanying datafile at hand. I couldn't find a formal specification for those index files either.

Perhaps one can recommend a library that can parse those index files, in python, C, or golang? It's hard to accept that there isn't a reusable block-gzipped index parser somewhere, other than the one buried deep in samtools/htslib/gatk code.

edit

I'd initially given a small motivating example making use of file offsets to access a subportion of file on a remote server. Commenters were quick to point out that samtools supported a series of protocols for doing exactly that: ftp, http, and s3 (which requires https, which requires that samtools be compiled with libcurl support) Thanks! That example was detracting, so I removed it.

samtools view does fetch a subportion of a file (potentially from a remote location), but it does a few extra things, like decompressing the chunks retrieved. With this question, I'm interested in manipulating the gz chunks directly.

So if the purpose of obtaining file offsets is to access subportions of remote files via protocols that samtools supports, as mentioned in this comment, then you can make your life easier and use samtools/tabix directly on the cli, prefixing your datafile with the corresponding scheme.

One side effect of using the CLI however, is that the blocks retrieved will be decompressed (e.g. samtools view). You'd have to re-compress them if what you're trying to do is create a smaller "bam" file.

Keep in mind that if you retrieve all zipped chunks which cover your range of interest, then you might be extracting more data than you asked for. The first and last chunk may contain records outside your region, just because their neighboring records happened to be compressed in the same chunk. You'd have to be careful to filter those downstream.

Specification for the tabix format

If you absolutely need a custom solution, the tabix data structure format, authored by Heng Li, is on the samtools website:

https://samtools.github.io/hts-specs/tabix.pdf

A more in depth discussion of how the data structure works can be found here:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3042176/

The tree-shaped binning index is a bit hard to understand from the text alone, but the tabix paper mentions that their approach is almost identical to the binning algorithm in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC186604/pdf/X6.pdf (page 1003). That paper has a nice illustration which makes things clearer.

Operation

Note: This is from my understanding of how the indices interoperate.

You don't need the data file to retrieve chunk offsets, you can do the work with just the tabix file. But it's complicated, and probably error prone if you're implementing this yourself.

Let's assume the query interval is Chr01:100-Chr02:2000.

Each sequence (Chr01, Chr02, etc) has its own pair of indices (binning index + linear index) in the tabix file. So you would break the query into two queries: Chr01:100-999999999999 and Chr02:0-2000, and concatenate the two results.

Within one sequence, to find all the blocks for open interval [X,Y), you would use the binning index to retrieve all the bins overlapping this interval (There is C code to do this in the data structure pdf). Then, you need to go through all the possible bins selected, and filter chunks. Not all chunks in a given bin will be overlapping with X, Y.

To filter chunks within a bin, you need the help of the linear index. There is one entry in the linear index per 16kbp of sequence positions. Each entry of the linear index contains the virtual offset of the leftmost record overlapping that 16kb window. I don't know what the entry contains if there is no record overlapping that 16kb window, presumably 0.

You would find the linear entry in which position X falls (X/16Kb), take that offset, call it O_x. You find the linear entry corresponding to position Y, and take that offset, call it O_y. When you go through the selected bins, filter out any chunk whose offset is < O_x, or whose offset >= O_y.

Assuming your datafile is sorted, the chunk range you're interested in is contained between the minimum offset and the maximum offset selected.

The virtual offsets are 64bit unsigned integers. The higher 48 bits of the virtual chunk offset corresponds to the file offset of the datafile (the other 16 you don't need. They are the offset of the record within that decompressed chunk).

Software

The software I have found typically mimics closely the hts_lib API. And hts_lib doesn't export a function which loads an index without also requiring a datafile. As far as I know.


nextflow: Filter outputs of a process How could I filter outputs of a process in the input of the next process? Filtering works fine in channel, but if I try to filter outputs I got compilation error. I tried as follows:
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}
    process step1a {
    input:
    tuple val(sample_name), path(bfiles) from data
    output:
    tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)
    script:
    """
    plink --bfile "${sample_name}" --chr 1-22 --out "${prefix}.step1b" --make-bed
    """
    }
    process step1b {
    input:
    tuple val(bim), path(bims) from step1a_bim.filter{bim, files -> name =~/*bim/}
    tuple val(sample_name), path(bfiles) from step1a_out
    output:
    path("${prefix}.step1b.snplist.txt")
    tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
    script:
    """
    awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""

It gives me
Compilation error
-cause: Unexpected input: '{'
process step1b {

Step1a gives me chr1.step1b.bed|bim|fam, chr2.step1b.bed|bim|fam.......chr22.step1b.bed|bim|fam (22 plink files)
I need to access prefix of those outputs to run plink and other tools to generate chr1.step1c.bed|bim|fam ...chr22.step1c.bed|bim|fam (another 22 plink files). run awk to generate chr1.step1b.snplist.txt .... chr22.step1b.snplist.txt (22 txt files).
Any help to filter the inputs so that I could use only *bim files from step1a_bim? Also how to select prefix of files in the input of processes that are from step1a_out to run plink? Does it make sense now?
Best Reagrds
Zillur

Instead of outputting all files and filtering for only bim, output only bims:
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)

Do
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.bim") into step1a_bim

So you can skip your complicated filter:
process step1b {
    input:
    tuple val(bim), path(bimfile) from step1a_bim
    output:
    path("${prefix}.step1b.snplist.txt")
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into
step1b_results
    script:
    """
    awk '{ if ((\\\\\$5=="T" && \\\\\$6=="A")||(\\\\\$5=="A" && \\\\\$6=="T")||(\\\\\$5=="C" && \\\\\$6=="G")||(\\\\\$5=="G" && \\\\\$6=="C")) print \\\\\$2, "ambig" ; else print \\\\\$2 ;}' "$bimfile" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"

"""

```


Supplementary aligments in VAF This question has also been asked on Biostars
I have a doubt, are supplementary alignment usually considered when the variant allele frequency is calculated? Thanks a lot.
In my case I have some Illumina short data edited with Crisp-Cas9 and I want to calculate for a specific position all the frequency of the 4 nucleotides. I would like to understand if the supplementary reads aligned in that position must be taken into consideration.
I've used pysamstats to calculate and it consider them, but my doubt is more if it right to consider them, if it is normally done or not. We are doing it on human and we expect only a change of a single nucleotide variant in the region if the editing works.
I'm interested in what is done generally, if there is a "rule" to consider or not the supplementary alignment when the AF is calculated.
In my case I have some Illumina short data edited with Crisp-Cas9 and I want to calculate for a specific position all the frequency of the 4 nucleotides. I was trying to understand if the supplementary reads aligned in that position must be taken into consideration.
GATK HaplotypeCaller and Mutect2 take them into consideration? I'm not looking for a specific variant caller, but if I want to take one variant caller as "reference", I would probably look at them.
Thanks a lot.

It depends on the tool and the use case but generally they should be.
It's not as important with short reads as only a small portion of reads have supplementary mapping (unless you're looking at breast cancer post chromothripsis, for example). But imagine the case where you're working with long reads (as I do). I have samples with a mean read length of 60kb. The human reference is structurally different to any given sample meaning that the chance of any given long read having supplementary mapping is high. If supplementary mapping didn't count towards allele counts then many would be artificially low. This is my use case, others are different.
I would count them unless you have some reason to suspect that they are are artefact of the editing and or the mapping. Sorry to never give a straight answer but that's bioinformatics! From what you described I would guess 99% are real but there's normally some noise and if that's of concern your code should include some QC. I.e, consider if multiple reads support the sup mapping and their quality.


What are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped? I have indexed a gzipped reference with bwa: bwa index reference.fa.gz, which produces a series of other files reference.fa.gz.{amb,ann,bwt,pac,sa}. These are working fine with bwa alignment.

I have discovered that samtools does not take a gzipped reference, so I am planning to use an unzipped version of the reference for my workflow instead of dealing with two separate representations of the reference.

Are the files reference.fa.gz.{amb,ann,bwt,pac,sa} dependent upon the reference being gzipped? Do I need to reindex the unzipped reference.fa with bwa, or can I just rename the current files to remove the .gz portion of the filename?

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.


How to filter blast results from blast run against virulence factor database I  downloaded virulence factor database from here , to predict virulence genes from few genomes that my lab  sequenced. I  performed local blast on the database with default setting and found around 4000 genes. I want to filter this result. I am not sure how to proceed. There are evalue, pident (percent identity), qcovs (query coverage) etc that I can tweak and I did tweak. For example, if I set pident above 40% the qcovs automatically goes close to 100%. That gave me indication above 40% is pretty significant and biologically relevant. But I am not sure what percent I should consider for treating the query as having the same function as the subject -- in this case that function is virulence.

Need some guidance.

Edit1: I ran blastp

Edit2: I am blasting protein (cds predicted from genome by annotation tool prokka). My question here is what is the most reasonable way to filter putative virulence factor using blast.

The answer is it depends. 40% pident with ~100% qcovs (assuming length of protein is not very small) tells you that the virulence factor (VF) and your prokka predicted protein are likely homologs, they belong to the same (or closely related) protein family. For some research studies, you might annotate your prokka proteins as 'putative' VFs.

On the other hand, in a clinical context, even 90% pident might not be enough evidence to say these predictions are real VFs. A single amino acid change alters the protein and therefore potentially the function. You could investigate bioinformtiaclly how the changes impact the protein structure and domains etc. Ideally validate it experimentally, e.g. for a new putative E. coli toxin, clone the protein and express it in a commensal E. coli strain and test the toxicity of the new strain.

In practise, think about what question you're asking and who you want to convince. Try numerous thresholds and see how robust your results are to threshold changes (4,000 VFs seems like a lot for a 'few' genomes). You could also try the VirulenceFinder database. When annotating antibiotic resistance genes, if you use 80% pident and qcovs thresholds with a minimum protein length then most researchers will accept your predicted annotations as 'conservative', maybe the VF community has similar arbitrary but pragmatic standards.


De novo motif discovery in protein sequences I am trying to build a features matrix to be used for Random Forest based classification. I'd like to add, as features, short motifs which are common to all the protein sequences belonging to a specific gene family.

I tried to use MEME but I'd like to know if there are also other good tools (online or command line). Keep in mind I have more than 10.000 sequences and that I prefer short motifs (3/4 aa long).

I think MEME is a good tool for your purpose, but there are others as well. I can think of InterProScan for example (although I am not sure if it is really de novo). Here a summary of available tools.


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

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


separate genes with expression values I have a file like this with 2 columns. I want to separate the columns based on plus and minus signs into separate files. How do I do this with awk or sed. Kindly help.

What do you want to do with NA values?

The below scripts ignore them (h3ab74's answer will print them to positiveExpression)

awk '$2 ~ /^[0-9]+$/  {print $0}' expressionFile.txt > positiveExpression.txt
awk '$2 ~ /^-[0-9]+$/ {print $0}' expressionFile.txt > negativeExpression.txt

These do a regex match to respectively, positive and negative numbers. If there's anything else in field 2 the line won't be printed.


Remove/delete sequences by ID from multifasta I have a fasta file like this:

>Id1
ATCCTT
>Id2
ATTTTCCC
>Id3
TTTCCCCAAAA
>Id4
CCCTTTAAA

I want to delete sequences that have the following IDs.

Id2
Id3

The IDs are in a .txt file, and the text file will be used to match and delete those sequences.

My output should be a fasta file like this

>Id1
ATCCTT
>Id4
CCCTTTAAA

But I want it with awk and/or sed and/or bash (no python or perl).
How can I do it?

If you want to learn how to do things with command-line tools, you can linearize the FASTA with awk, pipe to grep to filter for items of interest named in patterns.txt, then pipe to tr to delinearize:

$ awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff patterns.txt - | tr "\t" "\n" > out.fa

This will work on multiline FASTA.


Do you use SNAP for short-read mapping? I am calling SNPs from WGS samples produced at my lab. I am currently using bwa-mem for mapping Illumina reads as it is recommended by GATK best practice. However, bwa is a bit slow. I heard from my colleague that SNAP is much faster than bwa. I tried it on a small set of reads and it is indeed faster. However, I am not sure how it works with downstream SNP callers, so here are my questions: have you used SNAP for short-read mapping? What is your experience? Does SNAP work well with SNP callers like GATK and freebayes? Thanks!

GATK best practices are explicably meant to consume BWA MEM generated BAM.  Whilst SNAP may be faster, the Broad will not have tested it for compatibility with GATK as such you can't guaranty using it won't have unexpected consequences.

As such you'd be better off using BWA MEM because I assume accurately called variation is always better than fast and incorrectly called variation.  The main issue you'll have is ensuring shorter split hits and mapping quality are reported in the same way as bwa MEM -M which GATK/Picard is expecting.  Ultimately however you'd be better off posting this question on the GATK forum.

It's also worth noting that the soon to be released GATK 4 will utilise bwaspark which can distribute it's alignment processes across Apache Spark for increase performance.  Consequently I can't see SNAP being adopted anytime soon.


What is a poly(T) primer anchor sequence? I am reading Tang et al. (2009).

In the method description it is stated that:

  A single cell is manually picked under a microscope and lysed. Then mRNAs are reverse-transcribed into cDNAs using a poly(T) primer with anchor sequence (UP1) and unused primers are digested.

What is an anchor sequence in this case?

Tang, F. et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods 6, 377–382 (2009).

As shown in figure 1 of that paper, the anchor sequence is additional sequence that is added to the end of the polyT primer that provides a known and non-target binding site for additional PCR reactions. I gave a presentation about a month ago about the process that Oxford Nanopore Technology use in their strand-switching cDNA kit, which may be useful for understanding the process used by Tang et al.. Here are my strand-switch/PCR images from that presentation:

PolyT primer addition, with mismatched anchor sequence

Strand switch primer addition, with overhanging anchor sequence

PCR amplification, using primers that are complementary to the anchor sequences

In the case of the Tang et al. paper, they use a polyA-tailing method for strand-switching, but the process is essentially the same.


Obtaining all protein sequences with a particular domain architecture from Pfam I want to get the alignment of chain A of 1kf6 (PDB ID) from the pfam database here. This protein chain has two main domains (FAD_binding_2 and    Succ_DH_flav_C). In pfam there is a link to one of these domains and after clicking in one of the mentioned domains in the table below the page, another page comes at the top of which there is a link to architectures.

For example, if we click on "FAD_binding_2" in the table, there are 220 architectures that have this domain and after clicking on that button (at the top of page), a page opens here. In this page 7471 sequences with the following architecture: FAD_binding_2, Succ_DH_flav_C, resembles most to my chain (because it has both FAD_binding_2 and Succ_DH_flav_C domains).

If we want to see all 7471 sequences with that architecture we can click on the show button under that architecture and finally we can see all of that 7471 architecture resembling our chain A sequence. Down the page we can see other architectures that all have FAD_binding_2, but some of them do not have both FAD_binding_2 and Succ_DH_flav_C domains (or some of them have other extra domains that we are not interested in) and so we do not need their sequences and their alignments.

When we want to get the alignment of full sequences at the left side of the page, we can see that we can just get the alignment of all 15696 sequences and it seems that there is no way to get the alignment of our favourite 7471 sequences. I would like to know if there is any way to get the alignment or at least the fasta file of our favourite 7471 sequences (which have both of our favourite domains not just one of them)?

An alternative way is using the Hmmer website, and do a hmmsearch with the Pfam accession, forex. (PF02910). Navigate to the domain tab and find the domain architecture that you are interested in and press view scores. That way, only the sequences you are interested in will be downloadable under the tab downloads.

Step-by-Step guide:

1- Go to http://www.ebi.ac.uk/Tools/hmmer/

2- Search with hmmsearch http://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch

3- Use Accession search and type in PF02910

4- Navigate to Domain

5- To the right side of
9742 sequences with domain architecture: FAD_binding_2, Succ_DH_flav_C

Press View scores.

6- It should show (Your results have been filtered) and Query matches (9742)

7- Navigate to Download and download the sequences in the format of interest.

8- Repeat for domain architectures of interest.

I hope you find this helpful.

Kind Regards,

Sara


Correct way to assemble reference transcriptome - what's --samples_file option in Trinity? I have 101 samples from 9 tissues of a nonmodel species, coming either from SRA or sequenced by my lab. I want to generate a reference transcriptome with Trinity for further analysis. In order to do this, I downsampled the reads in order to have 20 to 30M reads for each tissue, for a total of 270M PE reads, which I combined in two files, one for left reads and one for right ones.

Now,  I came across the --samples_file flag in Trinity and I'm not sure if I should use it instead of concatenating the transcriptomes. Any help on good practice on this?

--samples_file flag is for DE analysis, you can put info about biological replicates in that file. You say that you want to make a ref transcriptome, so you don't need a samples file for DE analysis.

Why do you down sample? Are you not afraid you will not have enough reads to have a complete transcriptome now?


conda doesn't install latest version of snakemake I would like to use the --default-resources parameter for profiles which is available in later versions of snakemake.

To install snakemake, I created a new conda environment:

conda create -n snakemake -c conda-forge -c bioconda snakemake

The currently installed version is listed as 5.14:

$ conda activate snakemake
(snakemake)$ conda list | grep snakemake
# packages in environment at /home/uvi/be/hde/.conda/envs/snakemake:
snakemake                 5.14.0                        0    bioconda
snakemake-minimal         5.14.0                     py_1    bioconda

However, command line tool says it is on version 5.1.4:

(snakemake)$ snakemake --version
5.1.4

This version does not know the new parameter --default-resources yet. So it gives an error when using it

  snakemake: error: unrecognized arguments: --default-resources partition = "thin-shared"

I already checked that I'm running the executable from inside the environment, rather than a different, external one:

(snakemake)$ which snakemake
~/.conda/envs/snakemake/bin/snakemake
(snakemake)$ $(which snakemake) --version
5.1.4

Python (whether installed with conda or not) will prefer packages installed in your home directory unless you tell it not to. This is a little known secret about using a virtualenv, actually, since internally those tell python to ignore everything in your home directory. Within conda, there's long been a debate about whether python (and also R) should continue looking for packages in the normal places or only within a conda environment. Since it would likely break many people's workflows, it's been decided to allow the normal language features of looking for packages not just within conda envs but also in the normal system and user locations. For python in particular you can disable this as follows:

PYTHONNOUSERSITE=True snakemake ...

Then only the packages installed within your conda environment will be used.


Compare multiple alignment results' aligned bases I have aligned a nanopore data set to a reference genome with graphmap, minimap2 and BLASR. The alignment results are stored in BAM files.

I would like to do some concordance assessment, looking at the number of base pairs that are mapped to the same position by different aligner.

Are there any tools or ways that allow me to do this effectively?

I don't know of a single tool that can read multiple sets of BAM files and compare them. But you can use the Sam_Reader python3 module to create a csv file that summarizes each position in the genome for each of your 3 results and then compare them manually.

It should be installable with:

pip3 install sam_reader

And then generate the CSV with

from sam_reader import Sam_Reader
my_files = Sam_Reader('my_bam_file_or_folder/')
my_files.per_base_stats(write_file='new_file.csv')

*but FYI this module requires pysam, which does not work on windows


Arlequin files not able to converge beyond 2000 steps for some .arp files When I use .arp files in arlequin to run mismatch distribution it says that can’t converge beyond 2000 steps and also no graphs are produced.can someone suggests what can I do to overcome this problem

What you appear to have is a Bayesian MCMC for 2000 "steps" which isn't a lot and normally the term 'steps' is used for parsimony, so I'm a bit confused.

However, if a Bayesian analysis does not reach convergence, i.e. duplicate runs do not achieve the same results then the data space has not been solved so you get no result.

9 haplotypes does not sound much to me, it is possible there are too many sequences, it is also possible there is simply not enough variation (9 haplotypes ... how many alleles?).

Suggestions

Arlequin should provide a default convergence value which you can lower to force a result
Decrease the number of sequences but maintain the number of haplotypes
Try their example data set

Do you get any associated file such as posterior probability, any output file with a minus sign in front of it?


Read counts from BAM file I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:

samtools view -H IonXpress_input.bam

Output:

@HD VN:1.4  SO:coordinate   GO:none
@SQ SN:NR_039978    LN:984
@SQ SN:NM_130786    LN:1766
@SQ SN:NM_138932    LN:9293
@SQ SN:NM_033110    LN:2678
@SQ SN:NM_000014    LN:4678
@SQ SN:NM_144670    LN:5192
@SQ SN:NR_040112    LN:1201
@SQ SN:NM_017436    LN:2106
@SQ SN:NM_016161    LN:1771
@SQ SN:NM_015665    LN:1854

I guess these are from header. After these I see something like following:

@RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4

Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?

Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case).

For use in featureCounts my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).

If your data is properly aligned to hg19 your output should be something like this:

samtools view -H aligned2hg19.bam

@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10    LN:135534747
@SQ SN:chr11    LN:135006516
@SQ SN:chr12    LN:133851895
@SQ SN:chr13    LN:115169878
@SQ SN:chr14    LN:107349540
@SQ SN:chr15    LN:102531392
@SQ SN:chr16    LN:90354753
@SQ SN:chr17    LN:81195210
@SQ SN:chr18    LN:78077248
@SQ SN:chr19    LN:59128983
@SQ SN:chr20    LN:63025520
@SQ SN:chr21    LN:48129895
@SQ SN:chr22    LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566


Questions regarding the detection of variants (SNPs and Indels) in mammal Greetings,

I’m a newbie in bioinformatics working currently on my first project regarding detection of variants (SNPs and Indels) in mammal assembly genome. I used the following method to detect variants for non-model organism (Camelus dromedarius), and then predict variants effects:

1- Downloading two runs of raw sequencing reads in paired end mode using fastq-dump (version 2.9.3).

2- Raw reads were quality checked, then trimmed using Trimmomatic (version 0.36).

3- Clean reads were mapped to reference genome using BWA-MEM (version 0.7.15).

4- Merging SAM files to one big file, converting SAM to BAM, sorting reads by coordinate, and marking duplicate reads using Picard (version 2.20.1).

5- Create realignment targets using RealignTargetCreator, realign around those targets using IndelRealigner, calling variants for the Base Quality Scoring Recalibration BQSR process using HaplotypeCaller, apply BQSR using two steps of BaseRecalibrator then PrintReads, calling variants using HaplotypeCaller, previous processes were within GATK (version 3.8).

6- Calling variants again using BCFtools.

7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools.

8- Filtering variants using BCFtools.

9- Select variants SNPs and Indels in two separated files using GATK.

10- Predicting effects using SnpEff.

11- SnpEff database was built using reference genome + GTF file which was converted from reference genome GFF version 3 file using gffread.

Kindly, I would like to share my pipeline with those whom has an experience in related field, to review the above steps conducted on my task, whether it's accurate or still required additional steps, or rearrangement of the pipeline sequence?

In addition, I would like to share my results with you to make an assessment whether my results are suitable to be reported to my supervisor or not?

Your response is highly appreciated,

Kind regards,

This is close to what many people do. May I add:

between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this.
filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold.
For that purpose, some time ago I made this and this.