# Question
Title: Programmatically adding hydrogen and remove water to multiple PDB files
I'm trying to dock a ligand to several hundred PDB files (receptors). I thus need to prepare the those PDB files like removing water and adding hydrogen.

I can do that manually using Autodock or Pymol. But it's too troublesome. Is there a way to do it programmatically with Python or R or any command line procedure?

# Answer
> 2 votes
I believe PDB2PQR CLI will do the work wonderfully. Don't let the name trick you: PQR files are organized like PDB ones. Under the hood it runs propka, which is state-of-the-art for predicting a protein residues protonation state. The best part is that you can use PDB2PQR web server and they will give you the corresponding CLI arguments for each option you choose. It also automatically removes water (and other non-protein molecules).

You can also try looking into tutorials on how to setup proteins for dynamics on GROMACS, AMBER or OpenMM and stumble on their own tools for protonating stuff (respectively `pdb2gmx` (CLI), `protonate` (CLI), and `modeller.addHydrogen` (python API)), but as far as I know propka computes residues pKi in a knowledge-based manner and out-performs physics-based (if not just totally automatic, without guess) approaches like these.

For some context: acidic or basic amino-acids in proteins may have different pKi depending on their local environment (i.e. if they are in contact with solvent or not, and residues surrounding them), which is why "adding hydrogen to proteins crystallographic structures at a constant pH" is not as easy as it sounds.

---
Tags: python, protein-structure, docking, pymol, structural-biology
---


# Question
Title: how to get normcounts for singlecellexperiment object?
I need to perform differential expression analysis using the scDD package from R, but I am not able to since I miss the **normcounts** assay in my SCE object (of course in the example they show, the assay is already there).

As I found out in the preprocessing steps, one can use this function to generate the missing assay, except that it does not work:

```
# this is the SCE object
> cTSubSet
class: SingleCellExperiment
dim: 19767 149
metadata(2): Samples Samples
assays(2): counts logcounts
rownames(19767): Xkr4 Gm19938 ... CAAA01147332.1 AC149090.1
rowData names(4): ID Symbol Type chrLoc
colnames(149): AAAGGATGTTCACGAT-1 AACCATGTCTTTCGAT-1 ...
  TGTGGCGTCTCTTAAC-1 TTTGGTTAGACGCTCC-1
colData names(17): Sample Barcode ... slm_1 slm_1.5
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

# setting parameters as suggested in the vignette linked above
> paR <- list(alpha=0.01, mu0=0, s0=0.01, a0=0.01, b0=0.01)
# setting zero.thresh to zero since I already removed unwanted genes
> cTSubSet <- preprocess(cTSubSet, zero.thresh=0, scran_norm=TRUE)

# the not so useful error:
Performing scran Normalization
Error in t.default(Data) : argument is not a matrix

```

I was also looking for other possible ways of obtaining such **normcounts** assay and I found this description:

> normcounts: Normalized values on the same scale as the original counts. For example, counts divided by cell-specific size factors that are centred at unity.

Still, I can't figure out what this may mean. Is there a way of obtaining such **normcounts**?

# Answer
I do not know this package but it seems to call `scran` so why bothering and not call `scran` directly, see for an extended discussion + code:

https://bioconductor.org/books/release/OSCA/normalization.html#normalization-by-deconvolution
> 2 votes

---
Tags: r, rna-seq, scrnaseq
---


# Question
Title: Nextflow: How to initialize channel from CSV file but only after CSV file created?
Process A runs a python script which generates a CSV file in some absolute file path say "/mnt/x/blah.csv".

Process B needs to consume the CSV file line by line, after process A has finished.

The trouble is that the channel to input to B is initialized at the beginning before "/mnt/x/blah.csv" is created by A so it fails. This occurs even if I collect the result of process A as input to process B (because channel is initialized regardless)

```
nextflow.enable.dsl=2

def get_the_csv_records() {
    Channel
        .fromPath("/mnt/x/blah.csv")
        .splitCsv(header:true, quote: '\"')
        .map{ it.TheFieldIWant }
}

process process_a() {

    output:
    val("I'm done!")

    '''
    python.exe ./make_the_csv_file.py --path=/mnt/x/blah.csv
    '''
}

process process_b() {

    input:
    val(previous_step_is_done_dummy_variable)
    val(the_field)

    '''
    echo ${the_field}
    '''

}

workflow do_a_then_b{
    process_a()
    process_b(process_a.out.collect(), get_the_csv_records())
}

```

This fails immediately with "No such file: /mnt/x/blah.csv" because the channel is initialized before process\_a has even started.

Same behavior occurs when using concat:

```
    process_b(process_a.out.collect().concat(get_the_csv_records()))

```

Question:

1. How can I make process\_b() only read the input channel from the CSV file after it is created?
2. Alternatively, can I make process\_b() trigger when process\_a() is finished some other way?
3. Alternatively, can I add an intermediate step which reads process\_a.out.collect(), then reads /mnt/x/blah.csv and outputs it into some existing channel which process\_b() takes as input?
4. Another approach - what is the way to create a channel which outputs the results of splitCsv, when the file is not known at compile time?

# Answer
Assuming you are using latest version of Nextflow which support the DSL 2 language, and assuming your python process output the file blah.csv, then the easiest way to move forward will be:

```
process process_a() {
    output:
    path("blah.csv")

    '''
    python.exe ./make_the_csv_file.py --path=blah.csv
    '''
}

process process_b() {

    input:
    val(the_field)

    '''
    echo ${the_field}
    '''

}

workflow do_a_then_b{
    process_a() \
        | splitCsv \
        | map{ it -> it.TheFieldIWant } \
        | process_b
}

```

You can also refer to the nextflow google group.
> 1 votes

---
Tags: nextflow, workflow-management
---


# Question
Title: Downloading PSSM file from NCBI PSI-Blast results after the UI update
It used to be easy to download off a results page of a PSI-Blast from NCBI the PSSM (scorematrix). Now I cannot seem to find an `asn` file that is a scorematrix —`Download>ASN.1` is the output of the search. Whereas before it was in "format results" IIRC. I was after this to check something in the PSSM viewer or do further cycles of PSI-Blast with different settings —the latter especially handy given that at present the exclude taxa checkbox is troublesome. Where is the option gone?

# Answer
Have you tried clicking on the "Back to Traditional Results Page" link in the top right corner? If you're expecting a different interface, that might help work out where to find the right information.
> 1 votes

---
Tags: phylogenetics, psi-blast, pssm
---


# Question
Title: querying gdc metadata issues with resulting fields
I have been trying to access metadata from the Genomic Data Commons; however, not all of the information I want is showing up in the resulting tsv file. How do I get sample type and analyte type to show up in the tsv file? I have tried altering what I pass as a field (towards the top of the code), or adding a different filter, but nothing seems to change what shows up in the tsv file. Any help is appreciated!

```
import requests
import json

cases_endpt = "https://api.gdc.cancer.gov/cases"

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = [
    "sample_ids",
    "submitter_sample_ids",
    "files.cases.sample_ids",
    "files.cases.submitter_sample_ids",
    "files.cases.samples.sample_id",
    "files.cases.samples.sample_type",
    "files.cases.samples.sample_type_id",
    "files.cases.samples.submitter_id",
    "samples.sample_id",
    "samples.sample_type",
    "samples.sample_type_id",
    "tissue_source_site.bcr_id",
    "tissue_source_site.code",
    "tissue_source_site.name",
    "tissue_source_site.project",
    "tissue_source_site.tissue_source_site_id",
    "diagnoses.classification_of_tumor",
    "diagnoses.tumor_stage",
    "diagnoses.tumor_grade",
    "files.cases.diagnoses.tumor_grade",
    "files.cases.diagnoses.tumor_stage",
    "files.cases.diagnoses.classification_of_tumor",
    "files.cases.samples.tumor_code",
    "files.cases.samples.tumor_code_id",
    "files.cases.samples.tumor_descriptor",
    "files.cases.samples.portions.slides.percent_tumor_cells",
    "files.cases.samples.portions.slides.percent_tumor_nuclei",
    "samples.tumor_code",
    "samples.tumor_code_id",
    "samples.tumor_descriptor",
    "samples.portions.slides.percent_tumor_cells",
    "samples.portions.slides.percent_tumor_nuclei",
    "files.cases.samples.oct_embedded",
    "samples.oct_embedded",
    "files.cases.samples.intermediate_dimension",
    "files.cases.samples.longest_dimension",
    "files.cases.samples.shortest_dimension",
    "samples.intermediate_dimension",
    "samples.longest_dimension",
    "samples.shortest_dimension",
    "files.cases.samples.portions.is_ffpe",
    "samples.is_ffpe",
    "samples.portions.is_ffpe",
    "files.cases.samples.pathology_report_uuid",
    "samples.pathology_report_uuid",
    "exposures.weight",
    "files.cases.exposures.weight",
    "files.cases.samples.current_weight",
    "files.cases.samples.initial_weight",
    "files.cases.samples.portions.weight",
    "samples.current_weight",
    "samples.initial_weight",
    "samples.portions.weight",
    "files.cases.samples.composition",
    "samples.composition",
    "files.cases.samples.freezing_method",
    "files.cases.samples.time_between_clamping_and_freezing",
    "files.cases.samples.time_between_excision_and_freezing",
    "samples.freezing_method",
    "samples.time_between_clamping_and_freezing",
    "samples.time_between_excision_and_freezing",
    "files.cases.samples.days_to_sample_procurement",
    "samples.days_to_sample_procurement",
    "files.cases.samples.preservation_method",
    "samples.preservation_method",
    "samples.days_to_collection",
    "files.cases.samples.days_to_collection",
    "files.cases.samples.portions.analytes.analyte_type",
    "files.cases.samples.portions.analytes.analyte_type_id",
    "files.cases.samples.portions.analytes.aliquots.analyte_type",
    "files.cases.samples.portions.analytes.aliquots.analyte_type_id",
    "samples.portions.analytes.analyte_type",
    "samples.portions.analytes.analyte_type_id",
    "samples.portions.analytes.aliquots.analyte_type",
    "samples.portions.analytes.aliquots.analyte_type_id",
    ]

fields = ','.join(fields)

filters = {
    "op": "in",
    "content": {"field": "cases.project.project_id",
                "value": ["TCGA-BRCA"]},
}

params = {"filters": json.dumps(filters),
          "format": "TSV",
          "size": "12000"}

response = requests.get(cases_endpt, params=params)
content = str(response.content, "utf-8")

with open("addingsampletype5.tsv", "w") as file:
    file.write(content)

```

# Answer
> 0 votes
When I visit the `cases_endpt` page you linked directly on Firefox v88, I get a filterable JSON file. Looking for "type" only comes up with "disease\_type", not sample or analyte:

```
data
 hits
  0
   disease_type "Adenomas and Adenocarcinomas"
  1
   disease_type "Ductal and Lobular Neoplasms"
  2
   disease_type "Ductal and Lobular Neoplasms"
  3
   disease_type "Adenomas and Adenocarcinomas"
  4
   disease_type "Squamous Cell Neoplasms"
  5
   disease_type "Adenomas and Adenocarcinomas"
  6
   disease_type "Squamous Cell Neoplasms"
  7
   disease_type "Adenomas and Adenocarcinomas"
  8
   disease_type "Squamous Cell Neoplasms"
  9
   disease_type "Ductal and Lobular Neoplasms"

```

---
Tags: data-retrieval, api, python
---


# Question
Title: How can I create a combined barplot like this image?
I have data from different countries from two groups of samples (control and test). And, calculated three different indices alpha, beta and gamma. I have a large dataframe containing all the data. Here is a reproducible example of my dataframe (df):

```
df <- structure(list(
    SampleID = c("SRR490_profile", "SRR467268_profile",
                 "SRR467265_profile", "SRR467264_profile",
                 "SRR467263_profile", "SRR467258_profile",
                 "SRR467250_profile", "SRR467244_profile",
                 "SRR467242_profile", "SRR467227_profile",
                 "ERR710432_profile", "ERR710430_profile",
                 "ERR710429_profile", "SRR232646_profile",
                 "SRR232622_profile", "SRR232619_profile",
                 "SRR232617_profile", "SRR232616_profile",
                 "SRR232614_profile", "SRR232612_profile",
                 "HJ4386_profile", "HJ4382_profile",
                 "HJ4380_profile", "HJ4379_profile",
                 "HJ4377_profile", "HJ4375_profile",
                 "HJ4374_profile", "HJ4373_profile",
                 "HJ4372_profile", "HJ4370_profile",
                 "KLM3246730_profile", "KLM3246729_profile",
                 "KLM3246725_profile", "KLM3246722_profile",
                 "KLM3246721_profile", "KLM3246714_profile",
                 "KLM3246713_profile", "KLM3246712_profile",
                 "KLM3246711_profile", "KLM3246705_profile"),
    type = c("test", "control", "control", "control",
             "control", "test", "test", "test", "control",
             "control", "control", "test", "control",
             "control", "control", "control", "control",
             "test", "control", "test", "control", "control",
             "control", "control", "control", "test",
             "control", "test", "test", "control", "control",
             "test", "test", "test", "control", "test",
             "control", "test", "test", "test"),
    variable = c("gamma", "alpha", "beta", "gamma", "alpha", "beta",
                 "gamma", "alpha", "beta", "gamma", "alpha", "alpha",
                 "beta", "gamma", "alpha", "beta", "gamma", "alpha",
                 "beta", "gamma", "alpha", "alpha", "alpha", "beta",
                 "beta", "gamma", "alpha", "beta", "gamma", "alpha",
                 "beta", "gamma", "alpha", "beta", "gamma", "alpha",
                 "beta", "gamma", "alpha", "beta"),
    value = c(95, 90, 109, 112, 128, 109, 90, 87, 90, 99, 129, 101, 113,
              100, 111, 102, 99, 105, 85, 94, 133, 117, 116, 83, 113, 126,
              90, 128, 107, 121, 88, 109, 82, 42, 67, 100, 88, 101, 127,
              115),
    StudyID = c("swiss", "swiss", "swiss", "swiss", "swiss", "swiss",
                "swiss", "swiss", "swiss", "swiss", "australian",
                "australian", "australian", "australian", "australian",
                "australian", "australian", "australian", "australian",
                "australian", "thai", "thai", "thai", "thai", "thai",
                "thai", "thai", "thai", "thai", "thai", "english",
                "english", "english", "english", "english", "english",
                "english", "english", "english", "english")),
    row.names = c(NA, -40L),
    class = c("tbl_df", "tbl", "data.frame"))

```

Now, I want an image for each index containing barplots for two groups (control and test) and Y-axis showing the StudyID names (i.e., country names) like the following (this image shows plots for two indices alpha and beta):  Is this possible from my data? Can anyone give me a solution?

Thanks in advance

# Answer
> 2 votes
If you're not wedded to a dot and whiskers plot as the one show, I'd probably plot this dataset as a boxplot with ggplot:

```
ggplot(df,aes(x=value,y=StudyID,colour=type)) +
geom_boxplot() +
facet_wrap(~variable)

```

# Answer
> 2 votes
The ggplot2 function `stat_summary` or `stat_summary_bin` should get you most of the way there, something like this \[not tested\]:

```
library(tidyverse)
library(ggplot2)

df %>%
  ggplot() +
  aes(x=value, y=StudyID, colour=type) +
  stat_summary() +
  facet_wrap(~variable)

```

See `?stat_summary` for more information.

---
Tags: r
---


# Question
Title: Connection between Detected Genes and The Read Counts
I have been trying to understand the Seurat for analysing scRNA-seq data. It comes to my mind that the main data is organised in the Seurat object with rows as genes and columns as the cells, and the entries are the gene read count.

I am confused with the statement, "the connection between the detected genes in the cells and the read counts".

What is exactly the meaning of the detected genes in the cells biologically and its connection to the read counts? Can we say if the read count of gene A is 10 it means this gene is detected 10 times in a cell? How to interpret the detected genes biologically?

# Answer
*Answer from @ATpoint, converted from comment:*

A detected gene is one that has a count \> 0, and generally the deeper you sequence the more genes you detect until you reach a saturation which is being influenced by the number of biologically-expressed genes and the single-cell technology you're using, e.g. droplet vs plate-based and the library setup, e.g. end-tagged vs full length.

Counts cannot influence the nature of gene translation. A count is an observation; a translated gene sequence is a physical thing. Whether or not a gene can be translated into protein has no relationship to how often it is observed in existence.

Even with count zero (i.e. no detection in scRNA-seq), it does not mean that the gene is not present in the cell. NGS technology is a relative quantification, I doubt it is meaningful to try to **directly** make statements about the connection between read counts and biological meaning in the cell.
> 2 votes

---
Tags: scrnaseq, gene, seurat, reads, featurecounts
---


# Question
Title: How to make 3d model of a protein that not exist in PDB?
I'm working on nsSNPs(Non-synonymous single nucleotide polymorphisms) of protein AKAP13 whose PDB not available. I did one mutation in this protein and I want to show the effect by 3d modeling without using a PDB but don't know how to do. Can anyone help Please!!

# Answer
> 6 votes
## No PDB

Getting a crystal structure is hard work for crystallographers, even with a high throughput systems, so they make constructs with only known domains or regions of particular interest. What happened here is that only two parts have been cloned, expressed and crystallised. The rest may be crystallisable or solvable by cryoEM, but may be too disordered.

## Homology modelling

If the protein has a close homologue it can be threaded, that is it's sequences is "threaded" over the structure of the homologue. Two site provide these ready-made.

The former is very clear, has an NGL viewer, a handy circular map and the assessment is clear. For your protein it is showing four existing crystal structures, but these cover only a small part of it. The homology modelled section is small and not good here, but handy as the model has two numbers Qmean of –4 and identity of 23%, two red flags. So the next place to check are the following servers. They thread against different structures and hybridised. Missing parts are predicted ab initio (from the scratch). The latter parts look like spaghetti and cannot be trusted.

* Phyre2. \<1 day turnover. Do pay attention to the colour rails. If your region of interest in low confidence and is based on nothing, then do not trust it at all. Also there it is good to double check keep an eye on the sequence alignments of your position as misalignments can happen.
* I-Tasser. 3-5 day turnover. The CAPS server winner year after year. You 5 models representing each cluster. If they look very different and the metric is bad for the top, then you have very low confidence and must look at the all if even.
* EVFold. The amino acid co-ocurrence is rather interesting
* AlphaFold... Google keeps it secret. Nobody except them can use it.

Here is an example of a spaghetti model (not your protein):

The except for the red bar, which is modelled against homologues and forms a nice domain within the spaghetti ball at the bottom, all of the rest is *ab initio* and should not be trusted as the warning in the "Multi-template and ab initio information" section says.

A nice addendum to a model is Consurf, which maps conservation to the models as B-factors. In a crystal structure, the core is conserved and the surface is not. Metrics are nice, but human insight is better.

A big problem are protein that are part of a complex, because alone they are not structured and their structure is dependant on a neighbour. Models are great if they are threaded against a model that forms the same complex, eg. mouse to human, but not if they are hybridised.

In the case of multidomain protein that are large. It is best to split the protein up. Furtehrmore, I-Tasser and Phyre have length limits because the longer the protein to predict is, the higher the chance of error. One can determine where the domain boundaries may be —a Blast or a few iterations of Psi-blast will show conservation boundaries nicely. PFam has a lot to say about the protein architecture: enter link description here Regions of high disorder cannot be and should not be modelled as they do not form a structure free in solution —unless wrapped around another protein as mentioned.

## Wisdom to step away

However, the real question is not `how to` but `should one`. If one only get *ab initio* results, do not trust them at all —count your losses as opposed to embarking on a castle in the sky. For a SNV, check if it is close to a post-translational modification in PhosphoSitePlus or if there is a linear motif in ELM.

# Answer
> 2 votes
*Answer from @o-laprevote, converted from comment:*

Your protein includes intrinsically disordered regions (https://www.uniprot.org/uniprot/Q12802#family\_and\_domains). Showing a predicted structure may give a wrong idea of what it looks like: I'd consider simply showing the ordered part (of which a structure already exists) plus a "spaghetti" starting from its N-terminal residue.

In this case the ordered regions would be in the "interaction with ESR1" domain (1968-2338), which as you can see has been solved experimentally (PDB 4D0N), and the sequence spanning from 965-1431 (part of it has been solved, PDB 2DRN). What I mean by "spaghetti" is basically: render these two regions then draw curvy arabesques/lines between, before and after them to show it is disordered (starting from their N-terminal or C-terminal).

Such a cartoon is best done in a program that is shows it utterly differently —to not confuse the viewer as this is for display purposes only. Namely, render the parts in PyMOL, but then add the connecting lines in Photoshop or Powerpoint or MS-Paint/Preview.

---
Tags: python, 3d-structure
---


# Question
Title: Count of each sequence length from a fasta file with header using len function of python
I have a fasta file with 18 sequences inside. I want to extract the length of each sequence with header using the len function. My file looks like:

```
>DR-1
ATGTATGAGGGCAAGAAGACGAAGAACATGTTCCTGACCCGGGCCCTGGAGAAGATTCTG
GCCGACAAGGAAGTGAAGAAGGCGCATCACTCCCAGCTGCGCAAAGCGTGCGAGGTGGCC
CTAGAGGAAATAAAAGTTGAAACTGAGAAACAAAGCCCTCCTCATGGAGAAGCAAAAGCT
GGATCAGGCACCCTTCCACCAGTGAAATCAAAGACAAATTTTATTGAAGCAGACAAGTAC
TTTTTACCCTTTGAGTTGGCATGCCAGTCCAAATGTCCCCGCATAGTTAGTACATCTCTA
GATTGTTTGCAGAAACTTATTGCTTATGGGCACTTGACTGGCAGGGCTCCAGATAGTACG
ACACCAGGCAAAAAATTAATTGATAGAATTATTGAAACAATATGTGGCTGCTTCCAAGGT
CCTCAAACAGATGAAGGAGTTCAGTTGCAGATAATAAAGGCTTTACTTACCGCAGTGACA
TCACAACACATAGAAATCCATGAAGGAACTGTGCTGCAAGCTGTGAGGACATGCTACAAT
ATCTATTTAGCCAGCAAAAATCTCATCAATCAGACAACAGCCAAAGCTACTCTTACTCAA
ATGCTGAACGTTATCTTTGCACGCATGGAAAACCAAGCATTACAGGAGGCCAAACAAATG
GAAAGAGAAAGACATCGGCAGCAGCAGCATTTGTTACAGTCTCCAGTAAGCCATCATGAG
CCCGAGTCACCTCATCTTAGATATTTGCCACCTCAGACTGTTGATCATATAAACCAAGAG
CATGAAGGGGATCTTGAGCCCCAAACACATGATGTGGATAAAAGTCTCCAAGATGATACA
GAGCCTGAAAACGGATCTGATATTTCCAGTGCAGAAAATGAACAGACAGAAGCTGATCAG
GCAACTGCAGCTGAAACATTATCTAAAAATGATATCTTGTATGATGGAGACTATGAAGAA
AAGCCACTTGATATTGTCCAGAGCATTGTAGAAGAAATGGTGAACATCATTGTTGGAGAC
ATGGGAGAAGGGATGGCTATCAGTGCAAGTACAGAAGGCAACACTGGAACTGTAGAAGAT
GGTAGTGACAGTGAAAACATCCAAGCAAATGGAATCCCAGGAACACCCATTTCTGTTGCA
TATACACCATCGTTACCTGACGACAGGTTATCTGTCTCCTCCAATGATACTCAGGAATCT
GGAAATTCTTCAGGACCTTCACCTGGTGCTAAGTTTTCCCACATTTTACAAAAGGATGCC
>SK-1
TTTCTAGTGTTCAGGTCATTGTGTAAGCTGTCAATGAAACCACTGTCAGATGGACCACCA
GATCCAAAGTCTCATGAGTTACGGTCCAAGATTCTTTCATTGCAGTTGCTTCTGTCCATT
CTGCAGAATGCAGGACCTGTTTTCAGGACAAATGAGATGTTTATTAATGCTATCAAGCAG
TATCTATGTGTTGCACTCTCAAAAAACGGAGTCTCATCTGTTCCAGAGGTTTTTGAGCTT
TCTCTTTCCATATTTCTTACTTTGTTGTCAAACTTCAAGACACATCTGAAGATGCAAATT
GAGGTATTCTTTAAAGAAATCTTCCTGTACATTTTGGAAACGTCTACAAGCTCATTTGAC
CATAAGTGGATGGTTATTCAGACACTGACAAGGATCTGTGCAGATGCTCAGAGTGTGGTG
GACATTTATGTAAATTATGACTGTGACCTAAATGCAGCAAATATATTTGAAAGACTAGTA
AATGACCTGTCGAAAATTGCCCAAGGAAGGGGCAGCCAAGAACTTGGTATGAGTAATGTT
CAGGAATTGAGCCTGCGAAAAAAGGGGTTAGAATGTTTAGTTTCAATTTTGAAATGTATG
GTTGAATGGAGCAAGGATCAGTATGTGAACCCCAACTCTCAAACTACCCTTGGTCAGGAA
AAACCCTCAGAGCAAGAGATAAGTGAAGTTAAACACCCAGAGACTATAAACAGATATGGA
AGTTTAAACTCCTTAGAATCCACATCATCATCAGGAATTGGCAGCTACAGTACACAGATG
TCTGGCACTGATAATCCAGAACAATTTGAGGTCCTAAAGCAGCAAAAAGAAATAATAGAA
CAAGGGATTGATTTATTTAATAAGAAACCAAAGAGAGGAATACAATACCTCCAAGAACAA
GGGATGCTTGGCACCACACCCGAAGATATAGCTCAGTTCTTACATCAAGAAGAAAGATTA
GACTCTACTCAGGTGGGTGAGTTCCTGGGAGATAATGATAAATTCAATAAAGAAGTCATG
>SC-1
TATGCATATGTAGACCAACATGACTTTTCAGGAAAAGACTTTGTTTCAGCTCTTCGTCTA
TTTCTGGAAGGATTCCGGCTTCCAGGAGAAGCTCAGAAAATTGATCGATTAATGGAAAAA
TTTGCTGCAAGATACCTAGAGTGCAACCAAGGACAAACCCTTTTTGCTAGTGCAGATACT
GCTTATGTTTTGGCCTACTCAATTATCATGTTGACTACAGACCTTCACAGTCCACAGGTT
AAAAACAAAATGACAAAAGAACAGTACATTAAAATGAATCGTGGTATCAATGACAGTAAA
GATCTCCCTGAAGAGTATCTATCAGCCATTTATAACGAAATAGCTGGGAAAAAAATATCA
ATGAAAGAAACAAAAGAACTAACAATCCCTACAAAGTCAACTAAACAGAATGTAGCCAGT
GAAAAGCAAAGACGACTTCTGTATAACTTAGAAATGGAACAGATGGCCAAGACGGCTAAA
GCACTCATGGAAGCTGTGAGCCATGTTCAGGCACCTTTTACAAGTGCAACACACTTGGAG
CACGTGAGGCCGATGTTTAAGTTGGCTTGGACACCTTTTCTGGCTGCATTCAGTGTGGGT
CTCCAAGACTGTGATGATACTGAAGTAGCCTCTCTTTGCCTAGAAGGGATAAGATGTGCT
ATCAGGATTGCATGCATTTTCAGCATTCAGCTGGAAAGAGATGCATATGTCCAAGCACTA
GCCAGATTTACCTTACTTACAGTGAGTTCTGGTATTACTGAAATGAAACAGAAGAATATC
GACACAATAAAAACACTCATCACGGTGGCTCACACTGATGGAAATTATTTAGGAAATTCG
TGGCATGAGATCCTGAAGTGCATCAGTCAGTTAGAGCTGGCACAACTTATAGGAACTGGA
GTGAAGCCTCGGTACATTTCTGGGACTGTGAGAGGCAGAGAAGGGTCTCTCACTGGGACA
AAAGATCAGGCTCCAGATGAATTTGTTGGTCTAGGGCTAGTTGGAGGAAATGTGGACTGG
AAGCAGATAGCCAGTATTCAGGAGTCTATTGGAGAAACCAGCTCTCAAAGTGTGGTAGTT
GCTGTGGATAGAATATTCACCGGATCCACAAGGCTAGATGGAAATGCTATTGTTGATTTT
GTCCGTTGGCTATGTGCTGTGTCTATGGATGAATTGCTTTCTACTACACACCCAAGAATG
TTTAGTCTACAAAAAATAGTAGAAATATCATACTATAACATGGGAAGAATAAGACTGCAG
TGGTCTCGGATTTGGGAAGTTATTGGAGATCATTTTAATAAGGTTGGTTGCAATCCTAAT
GAAGATGTAGCCATATTTGCAGTAGACTCTCTGAGGCAGCTATCAATGAAGTTTTTAGAG
>SC-2
AAAGGAGAGCTTGCAAACTTCAGATTCCAAAAGGATTTCTTACGACCTTTTGAACATATA
ATGAAACGAAACAGGTCTCCAACAATTCGAGATATGGTTGTACGGTGTATAGCACAAATG
GTTAATTCTCAGGCTGCAAATATTCGTTCAGGATGGAAGAACATTTTCTCAGTATTCCAT
CTAGCTGCATCAGACCAAGATGAAAGCATAGTAGAACTTGCATTTCAGACAACAGGGCAC
ATTGTCACTCTTGTATTTGAAAAACACTTTCCGGCCACCATTGATTCTTTCCAGGATGCA
GTGAAGTGTTTGTCTGAATTCGCATGTAATGCAGCTTTCCCAGACACAAGTATGGAAGCA
ATCCGACTTATTCGCCATTGTGCAAAATATGTGTCTGATAGACCGCAGGCTTTCAAGGAG
TACACAAGTGATGACATGAGCGTGGCACCTGAAGACAGGGTGTGGGTGAGAGGGTGGTTC
CCGATCCTCTTTGAGTTGTCTTGTATCATCAATAGATGCAAATTAGATGTAAGAACCAGG
GGCTTAACAGTAATGTTTGAAATAATGAAAACATATGGCCACACTTATGAAAAACATTGG
TGGCAGGACTTATTCAGAATTGTTTTCAGAATCTTTGACAACATGAAATTGCCAGAACAG
CAGACAGAGAAAGCTGAATGGATGACAACAACTTGTAACCATGCCCTTTATGCAATATGT
GATGTATTCACCCAGTATTTAGAAGTTCTCAGCGATGTACTTTTGGATGATATTTTTGCC
CAGCTTTACTGGTGTGTACAACAAGACAATGAGCAGTTAGCACGATCTGGTACAAACTGC
TTAGAAAACGTGGTTATTCTGAATGGTGAAAAGTTTACTCTTGAAATCTGGGATAAAACG
TGTAACTGCACCCTGGATATCTTCAAAACCACAATTCCTCACGCGCTTCTGACTTGGCGT
>TA-1
CCCACCTCTGGCGAGGCCGAACCCCCATCTCCATCAGCTGTGAGTGAAAAGCCGTTGGAT
GCAATATCACAAAAATCTGTAGACATACATGATTCTATTCAACCAAGATCTTCAGATAAC
AGACAGCAAGCACCACTGGTTTCTGTTTCTACTGTGAGTGAAGAAGTTAGCAAAGTGAAA
TCCACAGCAAAATTCCCAGAACAAAAATTGTTTGCTGCCCTGTTGATTAAATGTGTTGTG
CAGCTGGAACTCATCCAGACCATCGACAACATTGTCTTCTTCCCAGCCACAAGCAAGAAG
GAGGATGCAGAGAACTTAGCTGCAGCCCAGAGAGATGCAGTAGACTTCGATGTTCGAGTT
GATACTCAAGACCAAGGAATGTACCGCTTTTTAACATCACAACAACTGTTCAAGCTACTG
GACTGCTTATTAGAGTCGCATAGATTTGCAAAAGCATTTAATTCCAACAATGAACAGAGG
ACTGCTCTGTGGAAAGCAGGCTTCAAAGGCAAGTCCAAGCCCAATCTTCTAAAACAGGAG
ACAAGCAGCCTGGCGTGTGGGCTGCGCATTCTCTTCCGGATGTACATGGACGAGAGCAGA
GTCAGTGCCTGGGAAGAAGTCCAGCAGAGGCTTTTAAACGTCTGCAGGGAAGCGCTGAGC
TACTTCCTCACTCTAACATCCGAAAGTCACCGGGAAGCCTGGACTAACTTACTGTTACTG
TTCCTCACTAAAGTTCTGAAGATCAGTGACAGCAGGTTTAAAGCTCACGCATCATTCTAC
TACCCTCTCTTATGTGAAATTATGCAATTTGACTTGATTCCGGAACTGCGTGCTGTTCTT
AGAAGATTTTTCCTGCGAATTGGAATAGTATTTCAAATATCACAACCACCTGAACAGGAA
CTTGGGATAAACAGGCAATGA
>DR-2
ATGCAGCCACCACCTCAGGCAGTCCCGTCTGGGGTGGCTGGACCACCCCCAGCTGGGAAT
CCTCGGAGCATGTTCTGGGCCAACAGCCCTTACAGGAAGCCAGCCAATAATGCACCAGTG
GCTCCCATAACACGCCCATTGCAGCCAGTAACGGATCCGTTTGCTTTTAATAGACAGACT
CTTCAAAACACACCAGTAGGAAGTTCATCCAAAAGCAGCCTGCCCAATTTGCCTGGCCCA
GCCCTCTCTGTTTTTTCTCAGTGGCCTGGTTTGCCTGTGACTCCCACAAATGCTGGGGAT
AGCTCCACAGGATTGCATGAGCCTCTGTCAGGAACTCTGTCACAGCCTAGAGCAGATGCC
AGTCTATTTCCTCCTGCATCTACCCCTTCTTCACTCCCTGGGCTGGAGGTGAGCAGGAAT
GCTGAGGCTGATCCCAGCTCAGGACATGAAGTCCAGATGCTGCCACATTCAGCTCACTAC
ATTCCAGGAGTGGGTCCTGAGCAGCCTCTCGGGGGCCAAATGAATGACAGTGGGTCCGGG
CCTGACCAACCCATGAATAGGCATGCCCCCCATGATGGTGCCGTGACCCACGCAGCATCT
CCTTTCCTCCCTCAGCCACAGATGCCAGGTCAGTGGGGGCCTGCACAAGGAGGCCCTCAG
CCCTCCTATCAGCATCACTCACCCTACCTGGAGGGACCTGTTCAGAACATGGGGCTCCAG
GCTGCCAGCTTGCCCCATTTTCCCCCTCCATCTAGTCTACACCAGGGACCTGGCCATGAG
TCTCATGCCCCACAAACATTTACTCCAGCATCCTTGGCCAGTGGTGAAGGGAATGAGATA
GTTCACCAGCAAAGTAAAAATCACCCATTGAGTAGCTTCCCCCCTAAACACACTTTTGAG
CAGAATTCCAGAATTGGGAATATGTGGGCAAGCCCAGAGTTGAAGCAGAATCCGGGAGTG
AATAAAGAGCATCTGCTAGACCCTGCCCATGTCAATCCCTTCACTCAGGGAAATAGTCCT
GAAAACCAGGCGCATCATCCCCCAGTGGCAGCAACTAACCATGCCCTGCAAGAAGCAGCC
TCTGGAGCTCTCTCAATGTTTTTTCAAGGGGAAGAGACAGAGAATGAGGAAAATCTCTCA
TCTGAAAAAGCAGGCCTCGATAAGAGATTGAACTTGGATAGCTTTTCTTCTACCTCCAGA
CTGGGCCATCCACCACCACCTGGAGCTAGTGGTGTCTACCAGGCCTTTCCTAGAGGCCCT
AGCAGTGAGGCTGCGCAGGAAGGAGATGCGCAACCTTATTTTTCGCAGTCTGTGGGCGTC
CGTCTTGACAAACAGAGCACTGTGCCCCCTGCTAATGATGCATGGGGTGATGTTCCGGGG
ACCGGGACCCGCTGTGCCAGTGGGCCACAGTGTGAGAATGTTGAGAACTTAGAATTTGTT
CAGAATCAAGAAGTTTTGCCCCGTGAGACCCTGAGTGTGGATCCCTTCCCCCTGAGTGAT
CAGATCCGATATGGGCCCCTTCCTGGGCCAGCTGCCTCCAGGCCTGCTACTGTGGGCCTC
ACCAGAGGTGGGGGGCTAAATCTTGAGGCCCCAGACACACCACTGCACCCTACACGACCT
GACAGTGTGTCATCCAGCTACAGCAGCCATAGCCACAGAAGTCCTCCTGGGTCAGCCAGA
CCCCAAGAATTGGTGGGCACATTCATTCAGCAAGAAGTTGGGAAACTTGAAGATGATACT
TCAGGGAGTTTTTTTAAGCAAATTGATTCTTCTCCTGTTGGAGGTGAAACAGATGAGGTG
ACTGGGAGCCAAAATTGCTGCAGCAGTCTGTCCCAGCCTTCAACCCCAAGCCCTCCAAAA
CCTACAGGAGTATTCCAAACAAGTGCAAACAGTTCTTTTGAGCCAGTAAAATCCCACTTG
GTTGGAGTGAAACCAGTTGAAGCAGATCGTGCCAACATGGTGGTTGAAGTGAGAGGGACC
CAATACTGTCCGAAGAAGCGCAGGGCAGCTGTTGCACCACCTGATGCTACCTCAGGAAAC
CTGGAACAGCCACCAGACAACATGGAGACCCCGTGTGCACCTCAGGCCTGCCCTTTGCCT
CTTAGCACCACTGGAGAGGCTGGGCAACTGGTTTCAAACACGGCAGGGACACCCCTGGAT
ACTGTGCGCCCAGTGCCTGATAAGAGGCCCTCAGCCAGGGCCCAGGGACCTGTGAAGTGT
GAGAGCCCAGCAACAACTCTGTGGGCACAAAATGAGCTGCCAGATTTTGGTGGCAATGTC
CTTCTAGCCCCTGCTGCCCCTGCACTTTATGTTCCTGTAAAACCAAAACCTTCAGAGGTT
GTTCATCATCCAGAAAAGGGAATGTCTGGACAGAAGGCCTGGAAACAAGGGTCTGTGCCA
>TA-2
CCTTTGCAGAACCAAGACCCTCCTGGTGCTTCTGAGAACCTTGAGAACCCTCCCAAGGTA
GGAGAAGAGGAGGCCCTTCCAGTGCAGGCAAGTTCTGGTTATGCGAGTTTGCTTTCTTCG
CCACCCACTGAATCTTTGCACAATCAACCTGTCCTGATTGCTCAGCCTGATCAAAGCTAT
AATTTGGCTCAGCCTATTAATTTTTCTGTGTCCTTACTGAATCCTAATGAGAAGAATCAG
TCCTGGGGAGATGCTGTGGTGGGGGAGAGATCCATAGTGAGCAACAACTGGGCTCTCGGT
GGTGATCCTGAAGAACGTGCTGCTCTGTCTGGAGTTCCAGCCAGTGCTGTCACTGGCGCG
TCTTTGCCTAGCAGTATCCCCCAGAATTGTGCCCCTCAAGGTTCGGGGTCTTCTGAAATG
ATTGCTAGCCAGTCTGCCAGTTGGTTGGTTCAACAGCTATCCCCTCAAACCCCTCAAAGT
CCGCATCCAAATGCAGAAAAGGGCCCTTCCGAGTTTGTTAGTAGCCCTGCTGGCAATACA
AGTGTGATGCTAGTTCCACCTGCAAGCAGTACTTTGGTTCCTAACAGTAATAAGGCAAAG
CACTCTAGTAACCAGGAAGAAGCTGTTGGCGCCCTAGACTTCACATTAAATAGAACTTTG
GAAAATCCTGTGAGAATGTACAGCCCATCTCCTTCAGATGGCCCAGCTTCTCAGCAGCCT
CTTCCCAATCATCCCAGACAATCTGGGCCTGGGCTGCATAACCAGGATCATTTCTACCAG
CAGGTAACAAAAGATGCTCAGGACCAGCATAGACTAGAAAGAGCCCAGCCAGAGCTGGTG
CCTCCTCGGCCACAGAATTCTCCCCAAGTGCCCCAAGCATCGTGTCCAGAGCCTTCAAAT
CCTGAAAGCCCACCCACACAAGGACAGTCTGAAAGCTTGGCCCAGCCACCAGCAAGTCCA
GCTTCAGTCAACACAGGCCAGCTGCTGCCTCAGCCACCACAGGCTTCCAGTGCATCCGTT
ACATCTACCAACTCCAGCCAGGCAGCTGTGCGGTCAGAGCAGCTGTGGCTACACCCACCA
CCTCCTAATACTTTTGGCCCAGCACCTCAAGACTTGGCATCCTACTACTATTACAGACCC
CTGTATGATGCCTACCAGTCTCAGTACCCTTCACCGTACCCGTCAGATCCTGGGACAGCC
TCGCTGTACTACCAGGATATGTATGGCCTGTATGAGCCTCGTTACAGACCCTATGACAGC
TCAGCATCTGCTTATGCCGAGAACCACCGCTATTCTGAGCCTGAGCGGCCCAGCTCCCGT
GCAAGTCACTATTCCGACCAGCTTGCTCCCAGGCAAGGATATCCTGAAGGATACTATAAT
TCCAAGAGTGGATGGAGCAGTCACAGTGATTACTATGCAAATTATTACTCTGGCCAGTAT
GATTATGGAGACCCCAGCCGCTGGGATCGTTACTATGGGTCTCGCCTTAGGGATCCTCGC
ACCTGGGACCGGAGGTACTGGTATGACTCTGAACATGACCCATACAGGAAGGACCACTAT
GCTTACAGTGACAGGCCTGAGAAATGTGATGATCACTGGAGGTATGACCCTCGCTTCACT
GGGAGCTTCGACGATGACGCTGAGATCCACAGGGACCCCTATGGAGAAGAAGCAGACAGA
CGCAGCATCCACAGTGAGCACTCGGCACGGAGCCTGCGCAGCACTCACAGTCTGCCCAGC
CGCCGCAGCAGCCTCAGCTCCCATTCACACCAGAGTCAGATTTACAGAAGCCACCATGTG
ACTGGCGGTTCCTTTGAGGCCCCTCATGCCCCAGGCTCATTTCATGGTGATTATGCCTAT
GGCACATATGCCAGCAATTTCAGCGGTGCCCATGGTTTTCCAGAGTACAGCTACCCTGCA
GACACCTCCTGGCCTGCTGTGGAGCAAGTTCCATCAAGACCGACCTCTCCTGAGAAGTTT
...
...
...
...

```

I want the output to look like :

```
>DR-1
1260
>SK-1
1020
>SC-1
1380
>SC-2
960
>TA-1
921
>DR-2
2400
TA-2
>1980
...
...
...
...

```

I am new to learning python for molecular analysis, so am having trouble with the logic and syntax behind it. I'm looking forward to a warm response from you. Thanks again.

# Answer
> 3 votes
Replace "string protein sequences.fa" with "your\_fasta\_file.fa" and your file should be at the same path that you created your code.

```
import Bio.SeqIO as IO
record_dict = IO.to_dict(IO.parse("string_protein_sequences.fa", "fasta"))
for key in record_dict.items():
    print(key[0],"\n ",len(key[1].seq))

```

Output:

```
9606.ENSP00000354718
  2817
9606.ENSP00000383118
  191
9606.ENSP00000275364
  381
9606.ENSP00000265563
  404
9606.ENSP00000265717
  418
9606.ENSP00000348461
  211
9606.ENSP00000249071
  192
9606.ENSP00000304283
  192
9606.ENSP00000400175
  193
9606.ENSP00000285735
  193
9606.ENSP00000339467
  191

```

I hope that's what you need, good luck

# Answer
> 0 votes
Don't forget to change the file name and write the valid path

```
with open("string_protein_sequences.fa") as file_content:
    liste=[]
    j=-1
    for seqs in file_content:
        if seqs.startswith('>'):
            j+=1
            print(seqs,lis[j])
            some=0
            i=0
        else:
            i+=1
            some+=len(seqs)
        liste.append(some-i)
    lis=[]
    for j in range(1,len(liste)):
        if liste[j]==0:
            lis.append(liste[j-1])
    lis.append(liste[-1])

```

Output:

```
 >9606.ENSP00000354718 (AKAP13)
  2817
 >9606.ENSP00000383118 (CDC42)
  191
 >9606.ENSP00000275364 (GNA12)
  381
 >9606.ENSP00000265563 (PRKAR2A)
  404
 >9606.ENSP00000265717 (PRKAR2B)
  418
 >9606.ENSP00000348461 (RAC1)
  211
 >9606.ENSP00000249071 (RAC2)
  192
 >9606.ENSP00000304283 (RAC3)
  192
 >9606.ENSP00000400175 (RHOA)
  193
 >9606.ENSP00000285735 (RHOC)
  193
 >9606.ENSP00000339467 (RHOG)
  191

```

# Answer
> 0 votes
you can use a dictionary as well:

```
fasta = {}
with open("test.fasta", "r") as file:
    for line in file:
        line = line.strip()
        if line[0] == ">":
            header = line
            fasta[header] = ""
        else:
            data = line
            fasta[header] += data

for k in fasta.keys():
    print(f"{k}\n{len(fasta[k])}")

```

---
Tags: python
---


title: alternative to HOMER's mergePeaks

text:
# Question
Title: alternative to HOMER's mergePeaks
Recently, our lab ran a 192 sample experiment through our ATACseq pipeline. In doing so, HOMER's `mergePeaks` told us that our 512 GB RAM server had too little memory to process all of the samples at once. So here we are trying to replace it with an alternative.

It seems that `bedTools merge` would be the way to go to replace it. However, going into this step, each sample has its own bed file of peaks called by `MACS2`. The `mergeBed` only takes in 1 input file that must be sorted by chromosome and peak start position before we run it. I can run a `cat` call on all of the samples followed by a `sort`, but those individual sample bed files are already sorted. Is there a way to run `mergeBed` with multiple input files at once? I tried `mergeBed -i 1.bed -i 2.bed > out.bed`, but that is equivalent to `mergeBed -i 1.bed > out.bed`.

In image format, I want the leftmost example, where the three peaks come from three separate files instead of all being described in the same file:

Additionally, bedtools does have a poorly documented function called `multiIntersectBed`that appears to run the `cat` and `sort` together, but is both slower and doesn't merge bookended sections--eg files look like this:

```
chr1 4658186 4658322
chr1 4658322 4658775
chr1 4658775 4658777
chr1 4671049 4671351

```

When they should look like this:

```
chr1 4658186 4658777 # This was 3 lines before
chr1 4671049 4671351

```

In summary, I have sorted bed files `1.bed` and `2.bed` from the same experiment, and I want to merge all bookended and overlapping peak sections between them. How do I do so in 1 step that doesn't involve HOMER's `mergePeaks`? Thank you!

EDIT: Here's the speed of my current solutions at n=192:

```
HOMER: Unstable

`multiIntersectBed | mergeBed`

real    3m53.230s
user    6m28.166s
sys     0m14.912s

`cat | sort | mergeBed`

real    0m24.860s
user    0m30.939s
sys     0m2.946s

`cat > tmp1; sort > tmp2; mergeBed > out; rm tmp?`

real    0m11.824s
user    0m20.710s
sys     0m2.021s

`sort -m | mergeBed`

real    0m17.905s
user    0m25.112s
sys     0m1.816s

`sort -m --batch-size=192 | mergeBed`

real    0m14.687s
user    0m25.589s
sys     0m0.747s

`sort -m > tmp ; mergeBed -i tmp ; rm tmp`

real    0m17.472s
user    0m15.251s
sys     0m2.212s

```

So... Followup questions, why does `sort -m` take longer than `sort | cat` especially when the sys time is indeed smaller? And why does piping take more time than creating temporary files?

# Answer
Via BEDOPS, you could combine `bedops --merge` with process substitutions sorting `N` files:

```
$ bedops --merge <(sort-bed 1.bed) <(sort-bed 2.bed) ... <(sort-bed N.bed) > answer.bed

```

The peak memory usage of a typical `bedops --merge` run should be around 10MB, if I remember correctly, as a result of using sorted inputs.

The memory usage of each `sort-bed` process can take up the size of each file. You can specify `--max-mem` and `--tmpdir` options to instead use disk space to sort files. You may want to do this if the sum of file sizes is equal to or greater than ~512GB, or however much system memory is available (typically less).

If your inputs are already sorted, or if your ATACseq pipeline is set up to create sorted files, then `sort-bed` does not need to be run again, and you can just use the merge operation on its own:

```
$ bedops --merge 1.bed 2.bed ... N.bed > answer.bed

```

The memory usage is then only what is needed for merging, as discussed above.

Using `bedops --merge` will merge adjacent ("bookended") intervals, so long as they are in sort-order, *e.g.*:

```
$ echo -e "chr1\t4658186\t4658322\nchr1\t4658322\t4658775" | bedops --merge -
chr1    4658186 4658775

```

Notes:

1. The file `answer.bed` will be sorted per `sort-bed` for downstream operations, such as `bedmap` association of merged regions to peaks. Additional sorting is not needed.
2. Using `sort-bed` will often be faster than Unix `sort` on BED files.
3. The limit of `N` on most platforms will be 1021 files. Most OSes have a limit of 1024 file descriptors, and three are reserved for standard IO. This soft limit can be dialed upwards on some hosts where the user has admin access.

References:

1. https://bedops.readthedocs.io/en/latest/content/summary.html#bedops
2. https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#merge-m-merge
3. https://bedops.readthedocs.io/en/latest/content/reference/file-management/sorting/sort-bed.html
> 1 votes

---
Tags: bedtools, homer
---


title: Pseudotime analysis on multiple samples

text:
# Question
Title: Pseudotime analysis on multiple samples
I am using `Monocle3` to run pseudo time analysis on samples at two different conditions (samples are paired). My question is, is it valid to combine all the cells from samples at specific condition when clustering, learning graph and ordering cells? or should I pick up only one sample per condition and run the analysis on?

My understanding is that pseudotime analysis uses the expression data of cells to try to infer the trajectory of cells, so I believe it might not be correct when I pool all the cells that might not have a true trajectory?

Thanks

# Answer
I would try it both ways (i.e. run pseudotime analysis on individual and combined samples). If you have the same transcriptional trajectories in multiple samples then you might be able to capture this via combining samples, and arguably the most interesting signatures are generalizable across samples. On the other hand, there might be batch effects that disrupt the signal of a combined analysis.

I think I understand your concern is that different samples have independent cells and so they should not form part of one trajectory. However, even within a sample, we are not capturing the true trajectory as this would require sampling the same individual cell multiple times at different time points, which is not currently possible because scRNA-seq requires lysing the cells (see also live-seq).
> 1 votes

---
Tags: scrnaseq, monocle
---


# Question
Title: A database for RefSeq protein accession IDs
Are there any databases or tools that contain a list of all RefSeq protein IDs for the interest organism/microorganism, e.g. get a specific microorganism name and return all RefSeq Protein IDs that have been identified?

# Answer
How about using RefSeq itself? Say, you are interested in all E. coli proteins, you can use the query `"Escherichia coli"[Organism] AND RefSeq[Filter]` here.

If you know the Taxonomy identifier for the specific organism of interest to you, it can also be used. For example, if you are interested in only "Escherichia coli str. K-12 substr. DH10B", you can use the taxonomy identifier 316385 (obtained from here) in the query `txid316385[Organism] AND RefSeq[Filter]`.

To extract a list of accessions and save them to a file on disk, use the "Send To" menu at the top right corner, choose "File" as destination followed by "Accession List" as the format and download the file. Alternatively, you can use E-utilities API or the EntrezDirect tool on the Unix command line to download these data programmatically.

To address the question about E. coli proteins shared between different species... There's no need to parse FTP files to get all the RefSeq proteins for a bacteria species. While it's true that some RefSeq WP proteins have taxids above species level, like Enterobacterales WP\_003861789.1 where the identical protein is found in genomes from multiple species of the order Enterobacterales, those proteins are still retrieved by a species-level organism search. NCBI implemented some extra handling for the RefSeq WP proteins to treat the MULTISPECIES WPs like they belong to all of the species with genomes on which they're annotated. So of the 56947 WPs assigned as "Bacteria", 7401 are found by `"Escherichia coli"[Organism]`, and the species-level search will give you the full set of proteins.

You could also use the `AutonomousProtein2Genomic.gz` file from the RefSeq FTP path which has a mapping of how all WPs map to specific genomic sequences, but it's a really big file and it would take a long time to grep for "Escherichia coli".
> 0 votes

---
Tags: database, proteomics, refseq
---


title: Hg38 annotation tracks retrieval

text:
# Question
Title: Hg38 annotation tracks retrieval
I want to retrieve an annotation object (as GRanges) containing different UCSC annotation tracks such as CpG Islands, TSS, TFBS and promoters for hg38 Human reference genome. But "UCSC Genes" was phased out for hg38. What could be the best resource for this ?

Thanks

RT

# Answer
The UCSC Data Integrator can combine multiple tracks (CpG Islands, TFBS, cCREs, etc.) into one output annotation file, and it works with any assembly you choose. You can select which output fields you want, then download the tab-delimited file.

From there, you can use R, Python, and/or GNU tools such as `awk` and `sed` to coerce the data into a format which is useful for your purposes. You can then load it into R as a data frame and use the columns as arguments in the GRanges constructor to create a GRanges object.
> 1 votes

---
Tags: sequence-annotation, ucsc
---


# Question
Title: Java error when trying to run UBCG jar (java.lang.ClassNotFoundException)
I wanna install UBCG.JAR in ubuntu 18.04 but I cannot run the java File UBCG.JAR. when I type the code `java -jar UBCG.JAR` an error message appears.

```
Exception in thread "main" java.lang.NoClassDefFoundError: org/json/simple/parser/ParseException
at java.base/java.lang.Class.forName0(Native Method)
at java.base/java.lang.Class.forName(Class.java:398)
at org.eclipse.jdt.internal.jarinjarloader.JarRsrcLoader.main(JarRsrcLoader.java:56)

```

Caused by: java.lang.ClassNotFoundException: org.json.simple.parser.ParseException at java.base/java.net.URLClassLoader.findClass(URLClassLoader.java:471) at java.base/java.lang.ClassLoader.loadClass(ClassLoader.java:588) at java.base/java.lang.ClassLoader.loadClass(ClassLoader.java:521) ... 3 more

```
$ java -version
Openjdk version "11.0.7" 2020-04-14

```

# Answer
> 0 votes
You are probably using the wrong java version it should be 8

try "java -version"

# Answer
> 0 votes
I solved this problem but UBCG.JAR does not work with java 11 . It is necessary to install openjdk-8-jre and switch the default java.

**Install openjdk-8-jre**

```
$ sudo apt-cache policy openjdk-8-jre
$ sudo apt install openjdk-8-jre=1.8.0_292

```

**Switch the default java**

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

```

**Checking the default java**

```
    $ java -version

```

openjdk version "1.8.0\_292" OpenJDK Runtime Environment (build 1.8.0\_292-8u292-b10-0ubuntu1~18.04-b10) OpenJDK 64-Bit Server VM (build 25.292-b10, mixed mode)

**Run the UBGC.jar**

```
$ java -jar UBCG.jar
      ###############################
      #         UBCG ver. 3.0       #
      ###############################

```

---
Tags: java
---


title: "variable.features.n" in SCTransform

text:
# Question
Title: "variable.features.n" in SCTransform
Whats the role of "variable.features.n" in SCTransform function of Seurat? I try setting "variable.features.n" at three different numbers (2000, 5000, 10000). In the downstream cells are clustering in different groups?

# Answer
`variable.features.n` sets the number of features (you can think of genes in the case of scRNA-seq) you would like to use for the downstream steps such as clustering. Basically, using a gene that is expressed more or less at similar levels across different cell types would not be informative in terms of differentiating (for example via clustering) these cell types, this parameter helps you to choose the "most informative" (most variable across the whole data set) features.
> 1 votes

---
Tags: seurat, sctransform
---


title: Making colnames by number and character

text:
# Question
Title: Making colnames by number and character
I have a very big matrix I want to put colnames as cell1, cell2, et,.. but I get error

```
> colnames(new.dat)
NULL
> View(new.dat)
> dim(new.dat)
[1]     34 263086
> colnames(new.dat)=make.names(c(1: 263086, "cell"))
Error in dimnames(x) <- dn :
  length of 'dimnames' [2] not equal to array extent
> colnames(new.dat) <- paste0("cell", seq_along(new.dat))
Error in dimnames(x) <- dn :
  length of 'dimnames' [2] not equal to array extent
>

```

any help?

# Answer
> 3 votes
You should really test your code:

```
make.names(c(1: 5, "cell"))
[1] "X1"   "X2"   "X3"   "X4"   "X5"   "cell"

```

This is obviously wrong.

Better:

```
colnames(new.dat) <- paste0("cell", 1:ncol(new.dat))

```

---
Tags: r
---


# Question
Title: How do I generate sequence logo of ATGC against each position in R?
I am new in R sorry for the naive question. I want to generate the seqlogo for each position using R. Can anyone guide me? Thank you so much in advance.

```
1  GATGGGGTGGGCAACTTAGAACTT
2  ATAACTGT
3  CAGAGGCGAACT
4  GGGCTAGGG
5  GGGGGGTGGCAGCACTAGC
6  GGGTAATATG
7  GTCGAGCGGT

```

Best

# Answer
> 3 votes
It's good practice to google before asking on online forums. Searching for "r seqlogo" pointed me to the seqLogo BioConductor package, which might be really helpful to you: https://bioconductor.org/packages/release/bioc/vignettes/seqLogo/inst/doc/seqLogo.html

You will need your sequences aligned though.

---
Tags: r, sequence
---


# Question
Title: Pulling a numbered chromosome range file given a gene location from a lookup table ideally from command line or R
have a folder with roughly 1000 vcf files which have divided the human genome into chunks, the folder looks like this:

```
- main_programme_aggregated_chr1_1_194789237.vcf.gz
....
- main_programme_aggregated_chrX_14235_8759493845.vcf.gz
etc

```

I will divide this file name into:

```
*_chr[1-22, X,Y]_A_B.*

```

I have a separate tab divided file with the names of genes and their genomic coordinates (all the genes in the human genome on build 38) - file1:

```
chr1 19853 90835 WASH7P ENSG0000022454
...
chr10 38732 390853 JHBF ENSG0000578382
..
chrX 9532 908032 LOTY ENSG00005847

```

etc

I would like to use the gene name from the 4th column of file1 to find the relevant chunked file it falls within. So I would need to take column 1 of file 1, match it to chr\[1-22, X or Y\] and then find the chunk were the range of columns2+3 fall between (A-B) in the file name.

I do not need to go into the files themselves, just use their names.

I have been just manually looking up the gene coordinates and the relevant chunk but am sure this can be automated.

An example would be that I want to get the chunk file that the gene WASH7P falls under, so the relevant file would be main\_programme\_aggregated\_chr1\_1\_194789237.vcf.gz.

Ideally I could do this from the command line and the relevant chunk file then be inputted as a variable into some downstream workflows I have.

Would anyone be able to help? Many thanks for your time

# Answer
> 2 votes
Here is a solution with `bash`:

```
#!/bin/bash

# 1 - create an array of coordinate files

# use nullglob in case there are no matching files,
# not that important for this toy case
shopt -s nullglob

# make a(n) array/list of vcf file in the current (./) directory
coordinate_files=(./*.vcf.gz)

# going through each line of "file1"
while read line; do

    # 2 - get the gene related info from "file1"
    # 2.1 - change line to an array (line split by whitespace)
    line_array=( $line )
    # 2.2 - get the 4th column, gene name and chr, the first column
    gene=${line_array[3]}
    chr=${line_array[0]}
    # 2.3 get coordinates
    start=${line_array[1]}
    end=${line_array[2]}

    # loop through each file in the array created at step 1
    for ((i=0; i<${#coordinate_files[@]}; i++)); do

        # parse the info hardcoded in the filename
        chr_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $4}')
        start_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $5}')
        end_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $6}')
        end_in_filename=${end_in_filename%.vcf.gz}

        if [ "$chr" == "$chr_in_filename" ]
        then
            if [ "$start" -ge "$start_in_filename" ] && [ "$end" -le "$end_in_filename" ]
            then
                echo $gene
                echo ${coordinate_files[$i]}
            fi
        fi
    done

done <file1

```

The input is the two example coordinate files + "file1" in the question and the output is:

```
WASH7P
chr1 19853 90835 WASH7P ENSG0000022454
./main_programme_aggregated_chr1_1_194789237.vcf.gz

```

# Answer
> 2 votes
Below is a Python solution:

```
import os

files = {}

for root, dirs, files in os.walk('/path/to/vcf/files'):
    for file in files:
        # Example: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz
        chrom = file.split('_')[3]
        start = int(file.split('_')[4])
        end = int(file.split('_')[5])
        files[os.path.abspath(file)] = {'chrom': chrom, 'start': start, 'end': end}
    break

genes = {}

with open(file1) as f:
    for line in f:
        fields = line.strip().split('\t')
        chrom = fields[0]
        start = int(fields[1])
        end = int(fields[2])
        gene = fields[3]
        genes[gene] = {'chrom': chrom, 'start': start, 'end': end}

mapping = {}

for gene in genes:
    for file in files:
        if (genes[gene]['chrom'] == files[file]['chrom'] and
            genes[gene]['start'] >= files[file]['start'] and
            genes[gene]['end'] <= files[file]['end']):
            mapping[gene] = file
            break

with open('mappings.txt', 'w') as f:
    for gene in mapping:
        f.write(gene + '\t' + file + '\n')

```

---
Tags: r, linux, shell, genomicranges
---


# Question
Title: Gromacs molecular dynamic simulation throwing error during minimization
I am trying to investigate a membrane proteins influence on a lipid bilayer using molecular dynamics. I am new to MD and trying to run my first simulation. I have built the simulation input files using CHARMM-GUI and now I am using gromacs to perform the simulation.

CHARMM-GUI provides these commands to get gromacs to perform the initial minimization

```
gmx_d grompp -f step6.0_minimization.mdp -o step6.0_minimization.tpr -c step5_input.gro -r step5_input.gro -p topol.top
gmx_d mdrun -v -deffnm step6.0_minimization

```

The first command runs fine. The second command starts the minimization, but at step 15 it freezes.

```
...

Step 15, time 0.015 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms -nan, max 0.098209 (between atoms 16596 and 16597)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
  15683  15684   39.2    0.1109   0.1121      0.1111
  15683  15685   38.4    0.1115   0.1107      0.1111
  16596  16597   72.1    0.1111   0.1220      0.1111
  16640  16642   35.3    0.1116   0.1120      0.1111
  16643  16644   51.9    0.1108   0.1136      0.1111
  16643  16645   47.2    0.1114   0.1079      0.1111

```

In step 15, the rms is -nan, so I'm guessing that is the problem, and I am assuming that is coming from the difference between the "current" and "previous" is greater than constraint length like on this line here

```
16596  16597   72.1    0.1111   0.1220      0.1111

```

The problem is I don't have any intuition on how to fix it? The pdb file is a homology model from phyre without any modification other than adding a single chain id to all the residues. I'm using a 1 to 1 POPC and DOPC in both membrane layers with a CHARMM36m force field.

The previous steps also throw LINCS warnings about bonds exceeding 30 degrees, however they produce a real value for rms.

I've counted on tutorials to get this far and would appreciate some suggestions or good resources I can use to continue my simulation.

Thanks!

# Answer
> 1 votes
*Answer from @user1271772, converted from comment:*

Try on this site: https://mattermodeling.stackexchange.com/questions/tagged/gromacs. There's an entire GROMACS tag there.

*Answer from @vishal-kumar-sahu, converted from comment:*

Try running/resuming using `-cpi` the simulation again, I got NaN error and rerun kept continued.

---
Tags: molecular-dynamics, gromacs, membrane-protein, lipid-bilayer
---


title: RNASeq read coverage in protein space?

text:
# Question
Title: RNASeq read coverage in protein space?
1. I used `samtools depth` to find the per base-pair read coverage over a number of isoform contigs from my Trinity assembly.
2. I have also conducted a multiple sequence alignment of those isoforms using Clustal Omega DNA aligner.

In my final step, I want to assess the per-residue read coverage of the consensus amino acid sequence to demonstrate an idea of confidence in that sequence.

Does anyone have any tips on how to go about doing this and whether there are existing tools for such things?

# Answer
*Answer from @cephbirk, converted from comment:*

What I ended up doing was to choose the read coverage for the lowest covered nucleotide in each codon as each residue's read coverage.
> 1 votes

---
Tags: sequence-alignment, read-mapping
---


# Question
Title: How to version the code and the data during the analysis?
I am currently looking for a system which will allow me to version both the code and the data in my research.

I think my way of analyzing data is not uncommon, and this will be useful for many people doing bioinformatics and aiming for the reproducibility.

Here are the requrements:

* Analysis is performed on multiple machines (local, cluster, server).
* All the code is transparently synchronized between the machines.
* Source code versioning.
* Generated data versioning.
* Support for large number of small generated files (\>10k). These also could be deleted.
* Support for large files (\>1Gb). At some point old generated files can permanently deleted. It would be insane to have *transparent* synchronization of those, but being able to synchronize them on demand would be nice.

So far I am using **git** \+ rsync/scp. But there are several downsides to it.

* Synchronization between multiple machines is a bit tedious, i.e. you have to git pull before you start working and git push after each update. I can live with that.
* You are not supposed to store large generated data files or large number of files inside your repository.
* Therefore I have to synchronize data files manually using rsync, which is error prone.

There is something called git annex. It seems really close to what I need. But:

* A bit more work than git, but that's ok.
* Unfortunately it seems it does not work well with the large number of files. Often I have more that 10k small files in my analysis. There are some tricks to improve indexing, but it doesn't solve the issue. What I need is one symlink representing the full contents of directory.

One potential solution is to use Dropbox or something similar (like syncthing) in combination with git. But the downside is there will be no connection between the source code version and the data version.

Is there any versioning system for the code and the data meeting the requirements you can recommend?

# Answer
> 2 votes
I'm a bit late to answer this question, but we've developed pretty much exactly the system you describe, based on Git and git-annex. Its called DataLad, free and open source software, and currently mostly known within the neurosciences.

DataLad's core data structure are DataLad datasets: Git repositories with an optional annex for large or private data. Within a dataset, code and data can be version controlled in conjunction. You can stick to Git and git-annex commands for handling the dataset, but DataLad also offers a core API that aims to harmonize and simplify version control operations. For example, in a dataset configured according to your preferences, a `datalad save` commits a file based on name (glob), file type, size, or location into either Git or git-annex (instead of `git add`, `git annex add`, and `git commit`). And likewise, a `datalad push` can push Git history and annexed data to your clones (instead of `git push` and `git-annex sync`). Everything is completely decentral, and DataLad integrates via git-annex special remote mechanism with a variety of hosting providers (S3, OSF, DropBox, ...), where you could store data if you like. A high level overview is here: handbook.datalad.org/r.html?about

> Unfortunately it seems it does not work well with the large number of files. Often I have more that 10k small files in my analysis. There are some tricks to improve indexing, but it doesn't solve the issue. What I need is one symlink representing the full contents of directory.

DataLad overcomes this issue by using Git's submodule mechanishm. Datasets can nested, and operated on recursively (to preserver a monorepo-like feel). With this mechanism, we've for example exposed the human connectome project dataset (80TB, 15 million files) to GitHub: github.com/datalad-datasets/human-connectome-project-openaccess

# Answer
> 15 votes
There is a couple of points to consider here, which I outline below. The goal here should be to find a workflow that is minimally intrusive on top of already using `git`.

As of yet, there is no ideal workflow that covers all use cases, but what I outline below is the closest I could come to it.

### Reproducibility is not just keeping all your data

You have got your raw data that you start your project with.

All other data in your project directory should never just "be there", but have some record of where it comes from. Data processing scripts are great for this, because they *already document how* you went from your raw to your analytical data, and then the files needed for your analyses.

And those scripts can be versioned, with an appropriate single entry point of processing (e.g. a `Makefile` that describes how to run your scripts).

This way, the state of all your project files is defined by the raw data, and the version of your processing scripts (and versions of external software, but that's a whole different kind of problem).

### What data/code should and should not be versioned

Just as you would not version generated code files, you should not want to version 10k intermediary data files that you produced when performing your analyses. The data that *should be* versioned is your *raw data* (at the start of your pipeline), not automatically generated files.

You might want to take snapshots of your project directory, but not keep every version of every file ever produced. This already cuts down your problem by a fair margin.

### Approach 1: Actual versioning of data

For your **raw or analytical data**, Git LFS (and alternatively Git Annex, that you already mention) is designed to solve exactly this problem: add tracking information of files in your Git tree, but do not store the content of those files in the repository (because otherwise it would add the size of a non-diffable file with every change you make).

For your **intermediate files**, you do the same as you would do with intermediate code files: add them to your `.gitignore` and do not version them.

This begs a couple of considerations:

* Git LFS is a paid service from Github (the free tier is limited to 1 GB of storage/bandwidth per month, which is very little), and it is more expensive than other comparable cloud storage solutions. You could consider paying for the storage at Github or running your own LFS server (there is a reference implementation, but I assume this would still be a substantial effort)
* Git Annex is free, but it replaces files by links and hence changes time stamps, which is a problem for e.g. GNU Make based workflows (major drawback for me). Also, fetching of files needs to be done manually or via a commit hook

### Approach 2: Versioning code only, syncing data

If your analytical data stays the same for most of your analyses, so the actual need to version it (as opposed to back up and document data provenance, which is essential) may be limited.

The key to get this this working is to put all **data files** in your `.gitignore` and ignore all your **code files** in `rsync`, with a script in your project root (extensions and directories are an example only):

```
#!/bin/bash
cd $(dirname $0)
rsync -auvr \
    --exclude "*.r" \
    --include "*.RData" \
    --exclude "dir with huge files that you don't need locally" \
    yourhost:/your/project/path/* .

```

The advantage here is that you don't need to remember the `rsync` command you are running. The script itself goes into version control.

This is especially useful if you do your heavy processing on a computing cluster but want to make plots from your result files on your local machine. I argue that you *generally don't need* bidirectional sync.

# Answer
> 10 votes
Your question is somewhat open, but I think it could prove an interesting discussion. I don't believe in many cases it is worth storing the data you have created in `git`. As you've noted, it isn't designed for large files (although we have `git-lfs`) and it's definitely not designed for binary formats such as `BAM`.

I'm of the opinion that `how` a file was created and `what` has been done to it since is key. Large files that took much effort to create should be mirrored somewhere (but not necessarily in a version control system). Other less-important (or less difficult to create) files that have been lost, clobbered or otherwise tainted can be regenerated as long as you know how they came to be.

For what it's worth, I've been working on a piece of software called `chitin` (self described as a shell for disorganised bioinformaticians). I also wrote a long blog post on why I thought this was a necessary project for me, but the main reason was that despite my attempts to organise my filesystem and make good archives of my experiments, over time I forget what my shorthand directory names meant, or exactly what program generated which data.

`chitin`'s goal is to automatically capture changes made to the file system during the execution of a command. It knows what commands to run to re-create a particular file, what commands have used that file and can tell you when and why that file was changed (and by who ;) too).

It's not finished (nothing ever is), but I feel that you might be going down the wrong road by wanting to store all your data and its versions when really, I think most people just want to know the commands that instigated changes. If the data history is important (and your code is well versioned), then you can simply check out any commit and execute your analysis to re-generate data.

# Answer
> 6 votes
First of all, kudos to you for taking versioning seriously. The fact that you're mindful of this issue is a good sign that you want to do responsible research!

For many bioinformatics projects, data files are so large that versioning the data directly with a tool like git is impractical. But your question is really getting at a couple of different issues.

* *How do I do my research reproducibly, and show full provenance for each data point and result I produce?*
* *How do I manage and synchronize my research work across multiple machines?*

**The short answer**:

* Archive the primary data.
* Place your workflow under version control.
* Version checksums of large data files.
* Use GitHub to synchronize your workflow between machines.

---

**The long answer**:

## Archive the primary data

As far as reproducibility is concerned, what matters most is the primary data: the raw, unprocessed data you collect from the instrument. If you are analyzing data that has been published by others, then write a script that will automate the task of downloading the data from its primary official source, and place that script under version control.

If you or a lab mate or a colleague produced the data and it is not yet published, then you should already have plans for submitting it to an archive. Indeed, most journals and funding agencies now require this prior to publication. I'd even go so far as to say the data should be submitted as soon as it's collected. Scientists worry a lot about having their data stolen and their ideas scooped, but statistically speaking, getting scooped is much less likely than nobody ever touching your data or reading your paper. But if you or an advisor insists, most data archives allow you to keep data private for an extended period of time until a supporting manuscript is published.

Putting (for example) Fastq files in a git repository is a bad idea for a lot of reasons. No hosting service will support files that big, git will be very slow with files that big, but most importantly git/GitHub is not archival! Use a proper data archive!

## Place your workflow under version control

Treat your raw data as read-only. Only process the raw data using scripts, and keep these scripts under version control. Vince Buffalo describes this well in his book Bioinformatics Data Skills. Check it out!

## Version checksums of large data files

If there are any data files that you want to track but are too big to place under version control, compute checksums and place *these* under version control. Checksums are very small alphanumeric strings that are, for all practical purposes, unique for each data file. So instead of putting that 5GB trimmed Fastq file or the 7GB BAM file under version control, compute their checksums and put the *checksums* under version control. The checksums won't tell you the contents of your files, but they can tell you when the file contents change.

This should give full disclosure and complete provenance for every data point in your analysis. The workflow has a scripts/command for downloading the primary data, scripts/commands for processing the data, and checksums that serve as a signature to validate intermediate and final output files. With this, anyone should be able to reproduce you analysis!

## Use GitHub to synchronize your workflow between machines

If your workflow is already under version control with git, it's trivial to push this to a hosting service like GitHub, GitLab, or BitBucket. Then it's just a matter of using `git push` and `git pull` to keep your code up-to-date on your various machines.

# Answer
> 5 votes
The Open Science Framework uses versioning for all files and is free to use: https://osf.io

You can integrate data or code from various sources such as github, dropbox, google drive, figshare or amazon cloud

You can also store files on their server using OSF data storage, but I do not know exactly what the file size limit is.

# Answer
> 4 votes
The way we deal with this is:

* All work is done in a single filesystem mounted on the cluster
* This file system is mounted on local machines via sshfs/samba (depending on the location of the current "local" machine on the network).
* Code is versioned with git hub
* All computation is carried out via light-weight automated pipelines. We use ruffus in combination with an in-house utility layer. The system doesn't really matter as long is it no more work to add another step to the pipeline than it would be to execute it manually.
* All questionably design decisions are encoded in configuration files. These configuration files, along with a very detailed log output by the pipeline (what was run, what was the git commit of the code run, what was the time stamp of the files it was run on, etc) and the initial input files are version controlled.
* The benefit of this is that code + configuration + time = output. It is not expected that the whole pipeline will be rerun everytime anything is changed, but the pipeline will tell you if something is out of date (it can use timestamps or file hashes), and it can all be run in one go before publication.
* Any other analysis is carried out in juptyer notebooks. These are version controlled.

To summarise, we don't synchronise because we only ever work from one disk location even if we use multiple CPU locations. We version control:

* Code
* Inputs, configuration, logs
* Juptyer notebooks

Log records the git commits used to produce the current outputs.

# Answer
> 3 votes
Using Git for version-controlling code is a good practice, but it does not lend itself well to versioning large data files. Manually syncing data across multiple nodes is asking for trouble, you want this syncing to either be handled automatically in a managed environment, or just keep the files on a single network-attached storage device.

One tool you might want to look into is Arvados, which is designed for syncing bioinformatics data and workflows across multiple machines. From the project website:

> Arvados is a platform for storing, organizing, processing, and sharing genomic and other big data. The platform is designed to make it easier for data scientists to develop analyses, developers to create genomic web applications and IT administers to manage large-scale compute and storage genomic resources. The platform is designed to run in the cloud or on your own hardware.

# Answer
> 3 votes
This answer will sort of only cover the big data parts, i.e. things \> 100MB ish, and only if your analysis pipeline ties in with the Python ecosystem. It will require a bit of learning

Try using quilt which is sort of like 'docker for data' (Github page).

```
$ pip install quilt
$ quilt install uciml/iris -x da2b6f5  #note the short hash
$ python
>>> from quilt.data.uciml import iris
# you've got data

```

Pros:

* There doesn't seem to be a filesize limit for public packages (haven't stress tested this though)
* Hashes your data to ensure reproducibility
* Support for versioning and tags

Cons:

* Storing private data starts from $7 per user / month up to 1TB of data
* Pretty much Python only at this point, with some community support for R

More information here.

---
Tags: reproducibility, data-management, git
---


# Question
Title: Voom function from limma package and Normalization on counts data
I know that Voom function from limma package from Bioconductor converts raw counts into log-CPM values and then Normalization is applied on that, with normalize.method argument.

I would like to know clearly for RNA-Seq data:

1. how this normalization is done (statistically) ?
2. normalization across samples and normalization across genes. But don't the know how it is done statistically.

# Answer
As you answered yourself @Death Metal, `voom` will by default not perform additional normalization.
However in virtually all cases you would want to do some kind of normalization at least to correct for differences in sequenced reads between the samples. This is why in the manual (page 71 right at the top) the `calcNormFactors` function from `edgeR` is used. This does a global `TMM` normalization for differences in library size.

Basically `voom` assumes that you did a normalization and will not do a second one because this could cause all kinds of problems.

In some cases you may prefer normalization other than `TMM`, for example if the number of genes detected in your samples is widely different. In this case is could be better to do `quantile` normalization that is available through the `normalizeBetweenArrays` function in `voom`.
However, this may over-correct your data and get rid of many interesting differences. In these cases I have found `qsmooth` quite helpful. Here quantile normalization is only applied across groups if the local distributions are similar enough.
> 2 votes

---
Tags: r, bioconductor, normalization
---


# Question
Title: What are some ways to error correct Oxford Nanopore long read sequencing?
I am sequencing long read genomic sequences to assemble MHC region haplotypes in a non model organism using Cas9 Sequencing Kit (SQK-CS9109) using flongle adaptors in the minION from Oxford Nanopore Technologies (ONT). However the error rate in genomic sequences results from ONT has been over 14%.

I know that there is error correction using Illumina sequencing and I was hoping anyone knows more about how this is accomplished?

Additionally, there are bioinformatic tools to error correct these errors? Has anyone accomplished diminishing this big error rate using bioinformatic tools?

# Answer
> 4 votes
How are you evaluating sequencing error rate? My most recent re-calls of 2017 sequences are demonstrating median single-read accuracy over 96%.

Before considering Illumina, it'd be worth it to do an initial correction using nanopore-only reads. This will make sure that the best results are obtained from the nanopore signal. First, make sure that the reads have been recalled using the latest highest-accuracy basecaller (currently guppy v5.0.11):

```
guppy_basecaller -c dna_r9.4.1_450bps_sup.cfg ... # for 9.4.1 reads
guppy_basecaller -c dna_r10.3_450bps_sup.cfg ... # for 10.3 reads

```

If you don't have any reference sequences available, reads can be corrected using the other reads by using canu:

```
canu -correct -p nonModel -d nonModel genomeSize=10m -nanopore reads.fastq

```

If you *do* have reference sequences available, then reads can be corrected using megalodon. I don't have any experience with megalodon, but have seen a couple of papers demonstrating that it has good accuracy, especially for methylation calling.

After all that, if you're still not happy with the accuracy, you can use Pilon for read correction using the Illumina reads. Pilon is used for correcting a *reference* sequence, so you'll need to create your haplotype assembly first before using Pilon. Even when using Pilon, I'd still recommend only using the Illumina reads for local INDEL correction. Pileup nanopore reads are very accurate for substitutions, especially when considering \[lack of\] strand bias when evaluating variants.

Update:

Given that you're interested in determining haplotypes across the MHC, which is something that performs poorly for short reads even for well-annotated genomes, I think you're better of sticking with reference-free correction, and categorising haplotypes while taking taking this into account. If you just need to be able to distinguish haplotypes, I recommend doing homopolymer compression on the resultant sequences, as hompolymer inaccuracies are the most common error mode for nanopore reads. Even if that doesn't get you all the way to perfect haplotypes, it might get close enough that something like CDHIT can be used for binning similar sequences.

HLA *typing* in humans can be done from nanopore reads (e.g. see here and here), but don't expect perfect base-level accuracy. What I'm suggesting is something similar to a typing process, but without an established reference, i.e. binning reads based on their similarity to each other.

---
Tags: nanopore, illumina, errors
---


title: Intrinsically disordered protein

text:
# Question
Title: Intrinsically disordered protein
*This question has also been asked on Biostars*

The position that I have to visualize is in the 2700 position and this position exists in two ranges in https://www.uniprot.org/uniprot/Q12802#family\_and\_domains the first range is Interaction with ESR1 and the second region is disordered, the first question is which range I should be based on? and the second is Why **intrinsically disordered protein(IDP)** can't be visualized in 3D?

# Answer
Asking how to visualise the disordered region of a protein is a bit like asking how to visualise the location of an electron.

An *intrinsically disordered protein*, by definition, has regions that are disordered. In other words, their location is variable under observation, so they cannot be precisely placed in a 3D model.
> 3 votes

---
Tags: 3d-structure, structural-biology
---


# Question
Title: Comparing the genetic similarity of human genome and mice/monkey genome - where can I find a list of genes?
I am trying to start a side project where I verify the claim that human and mice/monkey share 99% of their genes. I hope to learn the basics of doing genomic analysis.

I would like guidance about where I can find the sequences I need, and programming languages or applications that I should use. For the sequences, I would like for example, a list of all genes for each species. In the case that the programming language or application is foreign, some basic guidance would be helpful.

I have programming experience with matlab and python, and have learnt R previously, though forgotten. Would be able to learn a new language.

Alternatively, I think this may be a very large subject. If there is a link to a tutorial/lectures it will also be helpful.

Thanks in advance.

# Answer
For downloading lists of genes, together with associated features, I like using Ensembl Biomart:

http://ensembl.org/biomart/martview/

In this case, you can "CHOOSE DATABASE -\> Ensembl Genes", then "CHOOSE DATASET -\> Human Genes" or "CHOOSE DATASET -\> Chimpanzee genes" to get to a table selection.

Clicking on "Attributes" on the left hand sidebar brings you to a selector for different fields to add to the table.

Clicking on the "Results" button at the top left brings you to a form for choosing how results should be output.

Bear in mind that gene names are not likely to match. Finding the intersection of gene sets between species is a complicated process, usually requiring knowledge about gene copies, paralogs, orthologs, and sequence comparison.

It's also likely that the human genes are more annotated than the chimpanzee genes, and there might be genes present in the human list that are also in the chimpanzee *genome*, but have not yet been discovered.
> 1 votes

---
Tags: genome
---


# Question
Title: How do I get GO annotations for a list of UniProt IDs?
I have a list of UniProt ids that I want to get Gene Ontology annotations for. I need this information because I want this high-level information as an input to a neural network. The model I wish to develop is inspired by this paper: https://www.ijcai.org/Proceedings/2018/0468.pdf

Can anyone advise how I can programmatically access this information in the Gene Ontology Annotation (GOA) Database? https://www.ebi.ac.uk/GOA/

# Answer
> 4 votes
The link to the FTP for the GOA database files is listed on the GOA Downloads page. The file containing the mapping info you seek, `goa_uniprot.all`, comes in two formats, `.gaf` and `.gpa`. The README in the FTP directory linked provides details on the structure of the files.

Below is the first 20 lines of the `.gaf` file after uncompressing. Column 2 contains the UniProtID and column 5 contains the associated GO term. (Note that column 4 is empty in the snippet below)

```
>$ head goa_uniprot_all.gaf  -n 20
!gaf-version: 2.2
!
!The set of protein accessions included in this file is based on UniProt reference proteomes, which provide one protein per gene.
!They include the protein sequences annotated in Swiss-Prot or the longest TrEMBL transcript if there is no Swiss-Prot record.
!
!date-generated: 2021-02-12 21:14
!generated-by: UniProt
!go-version: http://purl.obolibrary.org/obo/go/releases/2021-02-08/extensions/go-plus.owl
!
UniProtKB   A0A009GUA7  J508_4179       GO:0003677  GO_REF:0000002  IEA InterPro:IPR006119|InterPro:IPR036162   F   Resolvase/invertase-type recombinase catalytic domain-containing protein    J508_4179   protein taxon:1310609   20201128    InterPro
UniProtKB   A0A009GUA7  J508_4179       GO:0006310  GO_REF:0000002  IEA InterPro:IPR006119|InterPro:IPR036162   P   Resolvase/invertase-type recombinase catalytic domain-containing protein    J508_4179   protein taxon:1310609   20201128    InterPro
UniProtKB   A0A009GUA7  J508_4179       GO:0000150  GO_REF:0000002  IEA InterPro:IPR006119|InterPro:IPR036162   F   Resolvase/invertase-type recombinase catalytic domain-containing protein    J508_4179   protein taxon:1310609   20201128    InterPro
UniProtKB   A0A009HCR2  J517_0313       GO:0005886  GO_REF:0000044  IEA UniProtKB-SubCell:SL-0039   C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0016020  GO_REF:0000043  IEA UniProtKB-KW:KW-0472    C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0016021  GO_REF:0000043  IEA UniProtKB-KW:KW-0812    C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0015031  GO_REF:0000043  IEA UniProtKB-KW:KW-0653    P   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0005886  GO_REF:0000043  IEA UniProtKB-KW:KW-1003    C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0055085  GO_REF:0000002  IEA InterPro:IPR003400  P   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    InterPro
UniProtKB   A0A009HCR2  J517_0313       GO:0022857  GO_REF:0000002  IEA InterPro:IPR003400  F   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    InterPro
UniProtKB   A0A009HLK9  J512_2607       GO:0006355  GO_REF:0000002  IEA InterPro:IPR000847|InterPro:IPR037402   P   HTH lysR-type domain-containing protein J512_2607   protein taxon:1310613   20201128    InterPro

```

---
Tags: proteins, uniprot, gene-ontology
---


# Question
Title: Help with DIA-Mass Spectrometry data analysis with several conditions (limma?)
I am trying to learn how to analyse normalised DIA-MS data and I am struggling with it ://

The original dataset I got is **(6 conditions (2 samples each)) with 3 technical replicates (total: 36 sample columns x 1128 proteins in rows)**. I removed the proteins with more than 80% of NA, and impute with kNN the rest of missing values.

Then I collapse the replicates doing the average of the 3 replicates (is this correct?)

**So now, I have a dataframe with 12 columns x 1072 rows.**

The question is... I am not sure how to perform the differential analysis from here (assuming I did well the previous steps lol)

Could someone help me with the tools or code to compare all the 6 conditions between them?

I read something about Limma package, but not sure if it can be used in this DIA-ms data and neither sure of how to use limma for it (I already check the vignette with no succeed)

UPDATE: I tried to use Limma as follows:

```
fit1 <- lmFit(norm.data_collapsed, design)
fit2 <- contrasts.fit(fit1, contrasts = contrast)
fit3 <- eBayes(fit2)
limma.res.01 <- topTable(fit3, coef = comb.contrast[1], n= Inf)

```

This returns me a non-sense values for logFC or significance:

Any help would be very welcome,

Best regards

# Answer
> 1 votes
Regarding the big FC values, I just realise the Limma lmfit function requires as input a log2 matrix... so I just transformed my data and now makes more sense. Still I have the doubt regarding the collapsing technical replicates and and if some one can verify the pipeline would be like that.

Thanks!!!

Fer

---
Tags: differential-expression, proteomics, limma, mass-spectrometry
---


# Question
Title: High percentage of poly A sequences in 10X chromium R2 read
I'm currently analyzing two samples of eosinophil cells isolated from mouse lung and the samples are of very different quality.
According to the Cell Ranger summary 56% of the reads can be mapped to the transcriptome in the first sample and only 32% in the second sample. Also the number of cells identified in Cell Ranger are very different (5000 for sample 1 and 600 for sample 2).

When I look at the raw data in FastQC I see that the second sample has a lot of poly-A sequences in the R2 read that obviously can't be mapped to the transcriptome.

Sample 1:

Sample 2:

Also in many cases there is the following sequence in front of the poly-A:
`AAGCAGTGGTATCAACGCAGAGTACATGGG`

To me it looks as if most of the reads do not contain an insert and are primer dimers. Is this connected to the lower cell count in the second sample?
Maybe some of you have already seen something like this and have a suggestion on how to avoid such a problem in the future.

Thanks a lot

# Answer
I actually figured out what was going on some time ago:

The over-represented sequence is the template switching oligo (TSO) that is used for the synthesis of the second strand (the example is for the v2 chemistry but the principle is the same for v3 and later).
This is a simplified illustration from the 10x manual:

As you can see the cDNA insert should be between the polyA capture probe and the TSO. In the second sample there are many cases with no insert and the TSO is followed by the poly A sequence. Actually the TSO should not be present in the final read after addition of Illumina adapters. But without any insert the fragmentation does not work and many reads still have this technical sequence.
So @merv was correct. Many reads contain no insert.

## Why did this happen?

Why were there many reads without insert in this experiment and why are some samples more heavily affected than others? As I mentioned, this were eosinophil granulozytes and these fellas contain many proteases, DNAses and RNAses to destroy invading pathogens. Unfortunately, they also destroy mRNA in single cell experiments.
Interestingly, the eosinophily in the second sample received an activating treatment and probably contained even more or more active RNAses.
> -1 votes

# Answer
*Answer from @merv, converted from comment:*

Perhaps it's an insert size issue. Can you get the post library construction QC info? You might also see this in the BAMs, with the first showing a full bell-shaped peaks, and the latter truncating the peaks at the cleavage site. (Most 10X datasets I've seen have very few reads traversing the 3' end cleavage site.)
> 1 votes

---
Tags: scrnaseq, ngs, 10x-genomics
---


# Question
Title: What is the origin of the 'source Gene ID' references given in the 'gene_presence_absence.csv' output of Roary?
I am learning to use Roary for preparing a pan genome for some lactobacillus strains. In the 'gene\_presence\_absence.csv' output of Roary (which I view in excel), a 'source Gene ID' is given for each of the strains included in the analysis. The Roary web page states that this is:

15. Presence and absence of genes in each sample, with the corresponding source Gene ID.

Is this 'source Gene ID' generated by Roary, or does it relate back directly to the .gff files I used as input? If I take one of these 'source Gene ID' references and `Grep` in the directory containing my .gff files I get no matches, so I am assuming it is generated by Roary...?

I had hoped that this reference would link back directly to the gene within the annotation file, but that seems not to be the case? Should/can I use this reference for anything?

# Answer
Looking at this again, I have now worked out (and should probably have noticed previously) that the 'source Gene ID' field from the 'gene\_presence\_absence.csv' file originates from the ID field of the .gff file for the sequence in question.

You can therefore easily tie the results from the 'gene\_presence\_absence.csv' file back to the annotation from the sequence in question, and therefore to the location of the nucleotide sequence for the gene in question.

So for instance, if you take a 'source Gene ID' from the gene\_presence\_absence file - say `GEICCBEJ_00950` \- and use it in a grep command referencing the genome from which it came, `grep "GEICCBEJ_00950" my_annotation_file.gff`, then you should see the line from the annotation corresponding to this sequence.

The 'source Gene ID' is also referenced in the pan\_genome\_reference.fa output from ROARY, which gives "single representative nucleotide sequence from each of the clusters in the pan genome (core and accessory)".
> 0 votes

---
Tags: genome, genomics, bacteria
---


# Question
Title: Impact of merging ChIP-seq runs of the same sample on PCR duplicates identification?
I'd like to a follow up question to this question related to merging fastq files for ChIP-seq.

Let's assume we have one sample library that is re-sequenced two or three times in order to achieve a desired read coverage. The resulting output of the sequencing is 2 or 3 fastq files for one individual sample.

If one has to mark duplicates (for example using Picard's MarkDuplicates) should the sub-samples be merged at the fastq level or at the bam file level (post alignment) after flagging duplicates before the merge?

If one takes the 2-3 sub-fastq files align them and then mark the duplicates in the bam files if the sequencing is shallow one could end up with not calling reads that are actually PCR duplicates of the library preparation. While in the other scenario, a fastq merge might result in reads to be marked as "duplicates" just because they appear 2-3 times, but maybe that read is a single reads that I sequenced 2-3 times, individually in each run.

My point is that I would expect this 2 different approached to lead to different outcomes. Probably this is a negligible effect in a overall analysis, but maybe it might actually have an impact in certain biological conditions (I'm thinking of heterochromatin silencing of repeat elements in the genome). I'm going to try it out as soon as I have a change, but I was wondering if this was discussed already somewhere else (I couldn't find much relevant) or if someone has some ideas on how to tackle this.

Let's assume that the reads are 50bp single end.

# Answer
If you are sequencing same library multiple times (as you state in the question), any read that appears in multiple sequencing runs is likely to have been generated via PCR duplication from the same original biological molecule, and therefore needs to be removed. Indeed, the chance of the read being a duplicate is the same if it was sequenced twice in the same lane, or it was sequenced in two different lanes, as long as it comes from the same library.

Thus, in this situation, you should either:

1. Merge the fastqs, map and MarkDuplicates or
2. Map the separate fastqs, merge the BAM files, and then MarkDuplicates.

You should definitely merge before doing either MarkDuplicates or peak calling.
> 1 votes

---
Tags: sam, fastq, chip-seq, merge, markduplicate
---


# Question
Title: Are variant calling files personally identifiable information?
This question is to be read in the context of data protection. It seems to be common understanding that the whole genome sequence of an individual is personally identifiable, non-anonymizable information: if you possess the whole genome sequence, you can without any doubt identify an individual. You also cannot anonymize it, because you would need to obstruct the genome sequence in some way, making it useless for any diagnostics or research.

But how about Variant Calling files (specifically referring to the VCF file format)? You reference against a known genome, but you neglect structural changes such as translocations, duplications and deletions. Can you still identify a person based on a Variant Calling file?

# Answer
> 9 votes
TLDR: yes! be careful with someone's genomic data!

There are two aspects to this question: can I find a match for a random VCF in a database of genomes (YES) or can I identify a subject who is not in a (public) database.

But even when it may not possible to identify the subject itself, using DTC genetic testing sites may enable you to identify (close or distant) relatives (e.g. read about the Golden State Killer). Identifying relatives brings you close to identifying the individual, and also those relatives may not be pleased with learning that the variants of someone with who they share (a small or larger part) of their genome are accessible by anyone.

We are further ahead for some traits than others, but for more and more human characteristics we will be able to predict e.g. eye color, height based on genome data, or at least estimate to some degree. Genetic ancestry and biological sex are trivial. Often genomic data is from a patient with a (rare) disorder, further narrowing down the options.

In addition, there is a difference between "can someone be identified today" (which is definitely possible but maybe tough) and "given that the variants are out there (online) forever can someone be identified within 10 years" and the answer to the latter part is a loud yes.

# Answer
> 8 votes
Yes, definitely identifiable. The combination of ~80 unlinked common SNPs can be fairly unique in the entire human population, let alone the whole VCF file.

EDIT: 30 in the original answer is an underestimate. We need ~80 unlinked common SNPs to uniquely identify every individual to a low false positive rate. Here is the derivation based on the approximate solution to the birthday paradox.

For simplicity, we assume each person has a haploid genome. Then the genotype at each SNP can only be 0 or 1. Suppose there are $n$ people in the world and we have $m$ sites. On the condition that $n^2\ll 2^{m+1}$, the probability of seeing any two people having the same binary genotypes across all $m$ sites is $$p(n,m)\approx 1-e^{-n^2/2^{m+1}}\approx n^2/2^{m+1}$$ There are $n\approx2^{33}$ people. Then $p(2^{33},72)\<1\\%$.

In reality, we are diploid, so each genotype has three states instead of two. This would lead to a smaller $m$ at the same false positive rate. On the other hand, selected SNPs are not totally random. This would lead to a larger $m$. It is safer to choose a larger $m$.

# Answer
> 3 votes
Not only are they theoretically identifiable based on the combination of SNPs as others have said (you are more or less identifiable from ~20 well chosen microsats) - but it's possible to be actually identified from your data using a long-range familial search, based on genetics and a bit of sleuthing.

https://science.sciencemag.org/content/362/6415/690.abstract

# Answer
> 1 votes
There should be someone at your workplace you can ask about this, but I would say, yes, since there can be SNPs in there that are private, or even a combination of public SNPs might be identifying, it should be treated as potentially identifying.

---
Tags: vcf, variant-calling
---


# Question
Title: Does it make sense to perform pathway enrichment analysis on a list of mutated genes from a VCF?
I have a list of mutated genes from a VCF. Does it make sense to perform a pathway enrichment analysis using PANTHER or STRING as examples to know which pathway is affected?

EDIT: to reply to @terdon

We sequence mostly human genome data using Illumina sequencers. As for the type of data, we perform mostly Whole-Exome Sequencing (WES), Whole-Genome Sequencing (WGS).

At the end, after performing all pre-processing, we obtain the VCF files from HaplotypeCaller using GATK. Then we proceed to annotate using VEP.

# Answer
> 1 votes
The concept of looking for enriched pathways makes sense, but you can't do it using the same tools you would use for gene expression data because the biases are different - longer genes are far more likely to contain variants by chance than short ones, and this needs to be accounted for in any enrichment analysis.

---
Tags: pathway, go-enrichment
---