# Question
Title: How drastic would the result be if I use fasta or reference assembly from ucsc and gtf from gencode?
There are difference annotation file for UCSC and gencode.

But if I use the reference assembly from UCSC and the GTF from Genocode or vice versa would my downstream results would be wrong?

# Answer
> 5 votes
**Never use genomes or annotations from UCSC**, they're poorly versioned and only recently with mouse and human have they even included all of the contigs. For fasta/GTF files from early in the GRCh38 release, you can tell whether you're using UCSC or Gencode by the presence/absence of `_random` contigs, which will only exist for UCSC. These were mostly later split into the actual contigs, so recent download from UCSC should more closely match what you find at Gencode/Ensembl. Further, that time predated UCSC beginning to adopt Gencode's vastly superior annotations, so if your GTF file has instances where the same gene ID is on either multiple strands or multiple chromosomes (this is obviously biologically impossible) then you have a UCSC GTF file.

In general, with early GRCh38 releases your only real issues will be with \_random contigs, which are a minority of the genome and don't have all that many genes. But really, you should be keeping track of the sources of your files and ensuring that they're compatible.

**Update**: I should expand a bit on my "Never use genomes or annotations from UCSC" comment. In point of fact the genomes themselves aren't so terrible. Early on UCSC had the bad habit of concatenating contigs together into `_random` "chromosomes", but they seem to have mostly kicked that habit as of late. Note, however, that there are no versions for their genomes. Since reference genomes continue to get updates over time (mostly through the addition of patches) the lack of actual versions means you have to manually check if a recently downloaded file matches what may have been downloaded either previously or by someone else. This has obvious consequences when it comes to reproducibility. The same issue occurs for annotations from UCSC, but they have the additional problem of historically having biologically incoherent concepts of genes. That is, they will contain the same gene in multiple places with multiple orientations, which will break many tools in both obvious and completely unclear ways. For example, DEXSeq will simply break with an error message if given a UCSC annotation, since they break biological plausibility. If you were to use these annotations files with deepTools, you wouldn't get an error message, but the resulting output would be only partial, due to the biologically impossible annotation effectively corrupting most obvious ways of storing annotation data in a data structure (i.e., you can no longer treat IDs as unique). This could have downstream ramifications on biological interpretation of results.

# Answer
> 1 votes
The main difference is in the way the chromosomes are named - UCSC uses the "chr" prefix (so chromosome 1 is "chr1") while in gencode the "chr" isn't used (so chromosome 1 is just "1"). Depending on your use case, this can obviously cause problems - if you're trying to match a locus (e.g. from gencode 1:1000002) between them, whatever tool you use is going to be looking in your aligned data for "1:1000002", but in your aligned data, it'll be named "chr1:1000002", so it won't match the two up.

---
Tags: rna-seq
---


# Question
Title: What is the `table` parameter in the UCSU ldHgGene tool?
I want to lift over the annotation from *E. coli* release 29 to *E. coli* release 42. I have followed the steps for an annotation liftover mentioned here. I got the chain file for the final step and ran:

`liftOver -gff Escherichia_coli.GCA_000967155.1.gtf Chromosome.chain conversions.gtf unmapped.gtf`

However, I get the following warning:

```
WARNING: -gff is not recommended.
Use 'ldHgGene -out=<file.gp>' and then 'liftOver -genePred <file.gp>'

```

I would like to heed the warning. However, I don't know how to use the `ldHgGene` tool. Here is its help page:

```
ldHgGene - load database with gene predictions from a gff file.
usage:
     ldHgGene database table file(s).gff
options:
     -bin         Add bin column (now the default)
     -nobin       don't add binning (you probably don't want this)
     -exon=type   Sets type field for exons to specific value
     -oldTable    Don't overwrite what's already in table
     -noncoding   Forces whole prediction to be UTR
     -gtf         input is GTF, stop codon is not in CDS
     -predTab     input is already in genePredTab format
     -requireCDS  discard genes that don't have CDS annotation
     -out=gpfile  write output, in genePred format, instead of loading
                  table. Database is ignored.
     -genePredExt create a extended genePred, including frame
                  information and gene name
     -impliedStopAfterCds - implied stop codon in GFF/GTF after CDS

```

From what I understand, I am trying to give the liftOver tool the annotation in GenePred format instead of GTF format. So I need to convert the GTF format into a GenePred format and that's what the ldHgGene tool is for. So the command should be something like:

```
ldHgGene Escherichia_coli.GCA_000967155.1.gtf -out=Escherichia_coli.GCA_000967155.1.gp

```

but this returns nothing. The `-out` option does mean that the `database` parameter is not needed but the tool still seems to need the parameter `table`.

What is the `table` parameter?

# Answer
I stumbled across this issue after getting the same warning message from `liftOver`. By using the `-out` option of `lgHgGene`, the `database` parameter is ignored but must still be provided. The `table` parameter should be the second column from your GTF file. Note that you must also specify the `-gtf` option (otherwise it expects GFF format). In my case, the GTF file looks like:

```
chr1    mm10_knownGene  exon    3073253 3074322 0.000000        +       .       gene_id "ENSMUST00000193812.1"; transcript_id "ENSMUST00000193812.1";
chr1    mm10_knownGene  exon    3102016 3102125 0.000000        +       .       gene_id "ENSMUST00000082908.1"; transcript_id "ENSMUST00000082908.1";
chr1    mm10_knownGene  exon    3206523 3207317 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3213439 3215632 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3205901 3207317 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
chr1    mm10_knownGene  exon    3213609 3216344 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
...

```

By default, `ldHgGene` outputs a leading 'bin' column which won't work with `liftOver`; this can be avoided with the `-nobin` option. So my command was:

```
ldHgGene -gtf -nobin -out=mm10.GENCODE_VM20.gp ignored mm10_knownGene mm10.GENCODE_VM20.gtf

```

The resulting genePred output can be fed directly into `liftOver -genePred`:

```
ENSMUST00000193812.1    chr1    +       3073252 3074322 3074322 3074322 1       3073252,        3074322,
ENSMUST00000082908.1    chr1    +       3102015 3102125 3102125 3102125 1       3102015,        3102125,
ENSMUST00000162897.1    chr1    -       3205900 3216344 3216344 3216344 2       3205900,3213608,        3207317,3216344,
ENSMUST00000159265.1    chr1    -       3206522 3215632 3215632 3215632 2       3206522,3213438,        3207317,3215632,

```
> 3 votes

---
Tags: format-conversion, sequence-annotation, ucsc, liftover
---


# Question
Title: Discordance of Gene Ontology data between Mouse Genome Informatics and Ensembl sources
In my first approach to the use of Gene Ontology, I have found a discordance if we consult either Mouse Genome Informatics (MGI) or Ensembl for these data.

Concretely, for the mouse genes Zic1 (NCBI 2271) and Zic2 (NCBI 2272), MGI assigns them up to 23 and 24 different GO terms, respectively.

However these numbers are different according to Ensembl data retrieved through R/Bioconductor. In particular, for Zic2 (Id 22772) the difference is very significant (24 terms in MGI vs. only 3 terms in Ensembl):

```
> my_ensembl = useMart ("ensembl", dataset = "mmusculus_gene_ensembl")
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22771', mart = my_ensembl)
        go_id
1  GO:0003676
2  GO:0046872
3  GO:0003677
4  GO:0005634
../..
22 GO:0008589
23 GO:0007389
24 GO:0021510
25 GO:0042472
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22772', mart = my_ensembl)
       go_id
1 GO:0003676
2 GO:0005634
3 GO:0016604
>

```

So, my questions are: Which is the reason for this discordance? Has it something to do with the different types of evidence to assign GO terms? Also, in order to proceed to analysis like GO enrichment, which of these two databases would be more appropriate?

# Answer
> 1 votes
Short answer, we get our annotations from GOA, which has the three listed. We're in touch with GOA to see why they differ from MGI.

---
Tags: bioconductor, ensembl, biomart, gene-ontology
---


# Question
Title: Get residue position from Pymol Selection
I can make a selection in Pymol using the resi parameter. Is there a way to get this resi data back out of a selection object? The only thing I found so far is how to get the position, as in the x/y/z coordinates from a selection.

For context, I have a table with a set of data points for each residue in the amino acid sequence of the displayed structure and I want to write a function where I can print the data based on the current selection. So my thought would be something like the following, where the content of the *'get\_residue\_positions\_for\_selection'* function would be what I'm asking about.

```
def show_selection_position_data():
    residue_data = {
        "header": ['variation', 'rmsf']
        , "1": ['high', 0.1]
        , ...
    }
    selections = cmd.get_names('selection')

    if sele not in selections:
        print("No selection")
        return

    residues = get_residue_positions_for_selection("sele")

    for residue in residue_selections:
        print(residue_data[residue])

    return

```

# Answer
> 2 votes
Ok, after some more searching I found the answer.

You have to use the **cmd.iterate** function to iterate over the default selection ("sele", first parameter). With the second parameter you can push the variables you need into an array for later access. In this case I wanted the residue number so I used the variable 'resv'. A full list of the variables exposed by the iterate function can be found on the pymolwiki

Note that you can use a function call as the second parameter and if you need a more complex selection than 'sele' you can use selector.process in the iterate function.

Here is my final working function. Note that I used a set for storing the residue ids, since the iterate function technically iterates over the atoms and so the same residue will likely appear multiple times in your selection.

```
def show_selection_position_data():
    residue_data = {
            "header": ['variation', 'rmsf']
            , 1: ['high', 0.1]
            , 2: ['high', 0.2]
            , 3: ['low', 0.1]
            , 4: ['high', 0.05]
            , 5: ['middle', 0.2]
            , 6: ['high', 0.1]
            , 7: ['high', 0.1]
        }
    selections = cmd.get_names('selections')

    if "sele" not in selections:
        print("No selection")
        return

    userSelection = "sele"

    stored.residues = set()
    cmd.iterate(selector.process(userSelection), 'stored.residues.add(resv)')

    for residue in stored.residues:
        if residue in residue_data:
            print(residue_data[residue])
        else:
            print("{}: no data".format(residue))

    return

```

---
Tags: pymol
---


# Question
Title: Why does DESeq2 convert numeric columns to factor during differential expression analysis?
I'm attempting to perform differential expression analysis using DESeq2 on a dataset of ~90 samples and ~20000 transcripts. I have a number of variables I'd like to test against, some of which are qualitative (factors) and some of which are quantitative (numeric). A sample of my data are below:

```
> head(data.phenotype.test)
# A tibble: 6 x 7
  Novogene_ID BORIS_ID     treatment queenness ovaries   elo sampletime
  <fct>       <fct>        <fct>         <dbl>   <dbl> <dbl> <fct>
1 W03BLGR     003W-blu-grn Qr12        0.0202    124.  1000  12
2 W03GR       003W-grn     Qr12        0.00210    66.8  923. 12
3 W03SL       003W-sil     Qr12        0.969     316.  1077. 12
4 W04BLOR     004W-blu-ora Qr12        0.00168    66.4  876. 12
5 W04GRSL     004W-grn-sil Qr12        0.175     210.   931. 12
6 W04YLWH     004W-yel-whi Qr12        0.984     302.  1192. 12

```

Now, let's say I want to test which transcripts are differentially expressed together with a numeric variable such as 'ovaries':

```
dds.transcript.ovaries = DESeqDataSetFromMatrix(countData = data.transcript.count.ovaries,
                                        colData = as.matrix(data.phenotype.ovaries),
                                        design = ~ ovaries)

```

This **should** work, but DESeq2 does something strange- it converts 'ovaries' to a factor, with the message *"some variables in design formula are characters, converting to factors"*. But there are no character values in the supplied data frame, and certainly not in the column 'ovaries'! This conversion to factors is a big problem for me because, having converted a numeric column to factors, DESeq2 unsurprisingly reports that differential expression analysis cannot be performed because there are an equal number of factor levels to samples.

I have found a temporary workaround: if I reduce the data frame to just the 'ovaries' column, DESeq2 no longer converts the numeric data to factor levels and I'm able to perform differential expression analysis as normal. The problem is that I'd eventually like to perform multivariate analyses, e.g. with a design like:

```
~ ovaries + elo + treatment

```

However, doing so isn't going to be possible if DESeq2 keeps converting the numeric columns to factors, which it persistently does unless the numeric columns are supplied alone (which of course they can't be for a multivariate design).

**What's going on here? Why does DESeq2 convert the numeric columns to factors, seemingly having misinterpreted those columns as characters?**

NB: I have replicated this issue with different data on a different machine and the result is the same, so I don't believe this an issue with my R setup.

Edit: Here's the structure of the df:

```
> str(data.phenotype.test)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   64 obs. of  7 variables:
 $ Novogene_ID: Factor w/ 64 levels "W03BLGR","W03GR",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ BORIS_ID   : Factor w/ 64 levels "003W-blu-grn",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ treatment  : Factor w/ 4 levels "Qc","Qr12","Qr3",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ queenness  : num  0.02017 0.0021 0.96853 0.00168 0.17529 ...
 $ ovaries    : num  124.1 66.8 316.1 66.4 209.9 ...
 $ elo        : num  1000 923 1077 876 931 ...
 $ sampletime : Factor w/ 3 levels "0","3","12": 3 3 3 3 3 3 3 3 3 3 ...
 - attr(*, "na.action")= 'omit' Named int  1 2 3 4 5 6 7 8 9 10 ...
  ..- attr(*, "names")= chr  "1" "2" "3" "4" ...

```

# Answer
Edited:

It's the "as.matrix" that is causing the problem. When I use that on my sample annotation data, and make dds with that, it turns my numbers into factors
> 1 votes

---
Tags: r, rna-seq, differential-expression, deseq2
---


# Question
Title: How to reverse complement the DNA sequences for given inverse/reverse coordinates?
I have the series of coordinates in `id.txt` file, whose coordinates sequences are in `genome.fasta` file. The coordinates of id.txt file is shown below,

```
Contig3:15-7
Contig2:5-10
Contig1:12-3

```

The genome.fasta file is shown below,

```
>Contig1
AAGGCCATCAAGGACGTGGATGAGGTCGTCAAG
>Contig2
ACGAAGCGCTCGCCAAGGCCGAAGAAGAAGGCC
>Contig3
GCTGCGGCGCTGATCCTGGCGGCCCGCGCCGAG

```

I have used the following codes to extract the sequences from `genome.fasta` based on the coordinates in `id.txt` by using the following command `xargs samtools faidx genome.fasta <id.txt > result.fasta` It has extracted only `Contig2:5-10`sequence. Because, the coordinates of `Contig2:5-10`is in proper order. But `Contig3:15-7` and `Contig1:12-3` coordinates are in inverse/reverse order, So samtools could not fetch those sequences.

I need to extract those sequences from `Contig1:12-3` and `Contig3:15-7` coordinates and also I want to `reverse complement` them. I have larger fasta file and plenty of coordinates sequences to be extracted, So I need to automate this process. So, please help me to automate this process.

The expected output is

```
>Contig3
GATCAGCGC
>Contig2
AGCGCT
>Contig1
CTTGATGGCC

```

Thank you in advance.

# Answer
> 1 votes
The `exonerate` tool comes with some very handy utilities for sequence manipulation. Here, we are interested in `fastasubseq` and `fasterevcomp`. First, install `exonerate` as described here. If you're using Debian or Ubuntu or one of their derivatives, it might be enough to run

```
sudo apt install exonerate

```

Then, follow these steps:

1. Create one fasta file for each sequence in `genome.fasta`

```
fastaexplode genome.fasta

```
2. A little preprocessing

```
$ awk -F[:-] '{
                if($2>$3){
                    start=$3-1;
                    len=$2-start" -"
                }
                else{
                    start=$2-1;
                    len=$3-start
                }
                print $1,start,len}' id.txt
Contig3 7 8 -
Contig2 5 5
Contig1 3 9 -

```

So, that will print the target contig's name, the start coordinate of the target sequence, the length of the target sequence and a `-` if it is on the reverse strand. The `fastasubseq` program I will use below starts counting from 0 and not from 1, so we also need to subtract 1 from the coordinates to get them right.
3. Now, feed that into `fastafetch` and `fastarevcomp`

```
tmpFile=$(mktemp);
    > subseqs.fa
    awk -F[:-] '{
                if($2>$3){
                    start=$3-1;
                    len=$2-start" -"
                }
                else{
                    start=$2-1;
                    len=$3-start
                }
                print $1,start,len}' id.txt |
    while read cont start len rev; do
        fastasubseq "$cont".fa $start $len > $tmpFile;
        if [[ -n $rev ]]; then
            fastarevcomp "$tmpFile" >> subseqs.fa;
        else
            cat "$tmpFile" >> subseqs.fa;
    fi
    done

```

This will produce the following file when run on your data:

```
$ cat subseqs.fa
>Contig3:subseq(6,9) [revcomp]
GATCAGCGC
>Contig2:subseq(4,6)
AGCGCT
>Contig1:subseq(2,10) [revcomp]
CTTGATGGCC

```
4. Finally, remove the per-contig fasta files

```
rm Contig*fa

```

# Answer
> 2 votes
Since you tagged python, here is a solution that solves this in a single script using Biopython:

```
import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start > end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')

```

Output on your example:

```
>Contig1
CTTGATGGCC
>Contig2
AGCGCT
>Contig3
GATCAGCGC

```

# Answer
> 1 votes
If you're okay with using a Perl script, you could convert your headers into stranded BED6 and use the following script `bed2faidxsta.pl` with `samtools`-indexed FASTA to get back stranded sequence:

https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640

This script looks at the sixth column of your BED file of chromosome, start, and stop positions, to decide whether to do a reverse complement of the queried sequence.

You could use `awk` to quickly build such a BED6 file:

```
$ awk -vFS='[:-]' '($0~/^>/){ if ($2 < $3) { strand = "+"; } else if ($2 > $3) { strand = "-"; } else { printf("Coordinates must differ\n"); exit 1; } printf("%s\t%d\t%d\t.\t0\t%c\n", $1, $2, $3, strand; }' in.fa > out.bed6

```

To set up indexed FASTA, here's a generic shell-based approach, which uses \[Kent utilities\]\[1\] `fetchChromSizes` to get a list of chromosomes for the build `mm9`:

```
$ mkdir mm9
$ cd mm9
$ while read chr; do echo "$chr"; wget -qO- "http://hgdownload.cse.ucsc.edu/goldenpath/mm9/chromosomes/$chr.fa.gz" | gunzip -c - > $chr.fa; samtools faidx $chr.fa; done < <(fetchChromSizes mm9 | cut -f1)

```

This can be adjusted for `mm10`, `hg19`, `hg38`, and so on.

Once you have indexed FASTA and the above Perl script, you could do something as easy as the following:

```
$ ./bed2faidxsta.pl --fastaDir="/path/to/indexed/mm9" < in.bed > out.fa

```

If you're working with a reference genome, the nice thing about this approach is that `awk`, Perl, and indexed lookups are fast. You can reuse indexed FASTA for other projects that also require sequence lookups by genomic interval.

---
Tags: python, fasta, samtools, shell
---


# Question
Title: Generating DNA sequences with constraints
I would like some advice on potential strategies to address the following problem.

I want to write a program that will generate DNA sequences that are optimized on two constraints based on an input genome:

* The GC content should be as close to the genome as possible.
* When taking all k-mers (of a predefined length) in the sequence, the sum of their frequencies in the genome should be as small as possible (i.e. it should be made of rare k-mers).

I currently score sequences using : $S = w\_{GC} * (1 - |GC\_{genome} - GC\_{seq}|) + w\_{Kmer} * (1 - \frac{mean(f(Kmer\_{seq})}{max(f(Kmer\_{genome}))}$

Where $w\_{GC}$ and $w\_{Kmer}$ are user-defined weights defining the relative importance of GC content and k-mer profile) and $f(Kmer\_{seq})$ is the frequency of k-mers in the sequence.

Of course I could just use bruteforce and generate many random sequences of appropriate GC content, score them all and pick the very best. But 1) that's no fun and 2) it's slow.

The first thing I attempted was to extend the sequence one basepair at a time. Each iteration has 4 candidates (with a new A, C, T or G added). I pick the candidate using their score ($S$) as probability weight. This "blind" method is too simple and does not yield satisfactory results.

I thought a dynamic approach could work well (filling a table and backtracking), but I am not sure this is applicable since the cost of extending the sequence would be dependent on all new k-mers generated and not be a fixed value (such as in sequence alignment).

I also thought something like reinforcement learning could be appropriate here, but I don't know much about it, so I am not sure.

I know all of this is probably overkill, but this is mostly for the challenge / learning :)

Which strategy would be most appropriate for this ? Did I miss something simple ? All advice / suggestions are welcome !

# Answer
You can use a $l$-order Markov chain. Here is the procedure:

1. Count $l$-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3.
2. Draw a $l$-mer randomly based on the distribution of $l$-mers. This generates the first $l$ bases.
3. Let $s$ be the last $(l-1)$ subsequence from the generated sequence. There will be at most four $l$-mers prefixed with $s$. It gives you a distribution of the next A/C/G/T base conditional on $s$. Randomly draw a base from this distribution. You get the next base.
4. Repeat step 3 until you generate a sequence long enough.

The generated sequence will have the same $l$-mer distribution to the original sequence. The GC content will be the same, too. In your case, if $k$ is small enough, you can let $l=k$. If $k$ is large, you have to choose a smaller $l$.

The procedure above doesn't consider the strand symmetry of DNA, but probably this doesn't matter much. It is also possible to adapt the algorithm to address this.

As to implementation, you can use python, but a high performance language like C will be tens of times faster at least.
> 4 votes

---
Tags: sequence-analysis, genomics, k-mer
---


# Question
Title: What is the purpose of 'Folder locations' in MaxQuant?
MaxQuant is a tool for proteomics quantification. In the settings under `Global parameters / Folder locations` are three folder locations that are possible to be changed. Thought the manual does not even mention these options. Is there any benefit in changing these locations? Or can it be bad to do that? For example, if I run multiple MaxQuant instances at the same time doing different searches?

What is the meaning of the `Fixed andromeda index folder`, `Temporary folder` and `Combined folder location` in MaxQuant? What is it used for and when would you (not) change it?

# Answer
These folders are usually generated in a default directory (where the raw files are). It can however be beneficial to set the paths to these folders manually.

> Temporary folder

* that is obviously where temporary files go, so having an SSD or other fast drive will benefit the I/O operations. Set that path to a fast storage partition.

> Combined folder

* all the output, including `txt` directory will go here. Setting it to fast storage might speed up writing the txt files. If you set it to the same path while running multiple instances of MQ it will obviously override output.

> Andromeda index

* this is where index files from Andromeda search are. By default these are inside a `combined/search` folder. Since these index files can get big, once generated by MQ they can be reused. Set this path if you want to reuse previously generated index files.
> 2 votes

---
Tags: proteomics, maxquant
---


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

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

```

My `din.fasta` alignment file shown below,

```
>org1
ATGCTA
>org2
TAGCAT
>org3
AGTAGT

```

I have tried the following code,

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

```

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

# Answer
> 4 votes
You were very close, you just need to supply the `offset` argument to `msaplot()`.

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

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

```

yields

Thank you for supplying a minimal and reproducible example!

# Answer
> 0 votes
I suspect the solution is to realign the two plots within "ggarrange", available here

The problem is stacking ggtree with geom\_tiplab. The ggarrange library is described on STHDA here

I have not used this library however, but it should be easy because you are working singly within R. I tend to use ggplot2 within Python via rpy2, which is much more difficult and restricted.

---
Tags: r, bioconductor, fasta, phylogenetics, phylogeny
---


# Question
Title: Find intervals that genes fall within their range
I am trying to run some frequency based stats to identify selective sweep in my study system using Rpackage "PopGenome".

To proceed with it I have split my genome data (whole genome sequencing) into 5MB chunks

```
cregions <- character(31)
cregions[1] <- "1-5000000"
cregions[2] <- "5000001-10000000"
cregions[3] <- "10000001-15000000"
cregions[4] <- "15000001-20000000"
cregions[5] <- "20000001-25000000"
cregions[6] <- "25000001-30000000"
cregions[7] <- "30000001-35000000"
....

```

I have a bunch of genes and I want to know each of my genes falls within the range of which intervals in the above list.

these are my genes:

```
chrom start stop gene_id
Chr01 107310787 107323633 gene2748
Chr01 113789595 113794851 gene2894
Chr01 159041389 159068535 gene4806
Chr02 2742629 2758060 gene4846

```

I'd like to get an out put like this

```
gene2748  window_2
gene2894  window_3
...

```

Any idea how I can do it in R?

# Answer
> 8 votes
The GenomicRanges package is a great place to start, something like this:

```
library(GenomicRanges)

# make a GenomicRanges set corresponding to the genes
genes.gr = makeGRangesFromDataFrame(df = genes, start.field = "start", end.field = "stop", seqnames.field = "chrom")

# turn your windows into a data frame and make a GenomicRanges set from that
windows = ...
windows.gr = makeGRangesFromDataFrame(df = windows ...

# find overlaps between the two sets
olaps = findOverlaps(genes.gr, windows.gr)

# output relevant fields
cat(paste(genes[queryHits(olaps),"gene_id"], windows[subjectHits(olaps),"window_id"]))

```

---
Tags: r, genomicranges, interval
---


# Question
Title: Error Code in WGCNA: Dealing with large data sets: block-wise network construction andmodule detection
I am comparing single block module detection with 2 block module detection in WGCNA.

Does anyone know what this relates to? I have rechecked and rechecked my codes for both the single and 2 module detection and cannot see an error. Thank you

I get the error code:

```
Error in plotOrderedColors(dendro$order, colors = colors, rowLabels = rowLabels,  :
  ERROR: length of colors vector not compatible with number of objects in 'order'.

```

Using the following code:

```
net= blockwiseModules(Th1_h5, power=6,
                      TOMType = "signed", minModuleSize=30,
                      reassignThreshold =0, mergeCutHeight = 0.25,
                      numericLabels = TRUE, pamRespectsDendro = FALSE,
                      saveTOMs = TRUE,
                      saveTOMFileBase = "PBCreponseTOMB",
                      verbose= 3
mergedColors = labels2colors(net$colors)

plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]

bwet = blockwiseModules(Th1_h5, maxBlockSize = 8000,
                        power=6, TOMType = "signed", minModuleSize = 30,
                        reassignThreshold = 0, mergeCutHeight = 0.25,
                        numericLabels = TRUE,
                        saveTOMs = TRUE,
                        saveTOMFileBase = "Th1_NR_R_TOM-blockwise",
                        verbose=3)
bwLabels = matchLabels(bwet$colors,moduleLabels)

bwModuleColors = labels2colors(bwLabels)

plotDendroAndColors(geneTree,
                    cbind(moduleColors, bwModuleColors),
                    c("Single block", "2 blocks"),
                    main = "Single block gene dendrogram and module colours",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang =  0.05)

```

# Answer
> 2 votes
The problem is that blockwiseModules split your data into (probably 3) blocks, because the default maxBlockSize is 5000 and smaller than the number of genes in your data. If you have enough RAM (16GB is recommended for 12k genes, but 8GB just might be enough), rerun blockswiseModules again with maxBlockSize raised to above the number of your genes. Since you were able to run the function with maxBlockSize=8000, hoepfully you'll be able to do the single block calculation.

---
Tags: r, wgcna
---


# Question
Title: What is the best way to start with Structural bioinformatics with zero prior experience?
My undergraduate thesis is on the topic "Ligand design from protein-protein interfaces" and my PI has asked me to learn the basics of Structural bioinformatics over the summer so that I'll be able to start with my project as I begin my next semester(last week of July). Can someone please suggest me top resources to learn the basics of structural bioinformatics. Though I have done projects in Genomics, and Systems biology, this will be my first time with structural bioinformatics, so, I'll really appreciate some good leads. Thanks!

# Answer
> 2 votes
As others have said, there is no short cut and you will need to spend some time with a good book.

I would say there are 3 broad areas of understanding you have to have as a structural bioinformatician:

* Solid understanding of the basic science. This includes molecular biology but also chemistry, which is often overlooked.
* Knowledge of the databases and data formats available in the field. You need to know where to go to get protein structures, protein sequences etc.
* Programming and algorithmic skills. You don't have to be a computer whiz but being effective in a high-level language (e.g. Python, R, Julia) and being able to parse data and carry out basic data analyses are essential.

If you want to take it further and have a bit more time then I would consider implementing something yourself such as a protein secondary structure predictor. Getting this working is challenging and would show you have competency in the above three areas.

# Answer
> 1 votes
A good starting point could be this book by Gu and Bourne:

Structural Bioinformatics, 2nd Edition

---
Tags: protein-structure, protein-protein-interaction, structural-biology
---


title: Occupancy of TFs with the target genes

text:
# Question
Title: Occupancy of TFs with the target genes
The occupancy of SMARCD3 in the target genes listed below. I want to see average, normalized ChIP-seq signal at the promoter proximal region (1000bp upstream and downstream of the TSS).

I have 4 different experimental conditions (overlayed in one plot, per target gene using different color codes) where we want to visualize the changes in occupancy of GATA1 under all the conditions at individual gene promoters. The following are our target genes

1. ABCB1
2. ABCC1
3. ABCC2
4. SNAI2

So how GATA1 occupancy is changed across all these target genes across 4 different treatment condition.

So far what i have done

1. Aligned the samples
2. Done peak calling
3. Annotated the peak

I did steps 2 and 3 using homer.

Now after all this how do I find the occupancy of a TF with its target genes? Do I get the information from my peaks annotated output or from the peak files?

# Answer
I would ignore peak calling for this and instead compute enrichment of ChIP/input for the genome (e.g., with deepTools or presumably homer) and then plot it for the genes of interest individually (e.g., using IGV or pyGenomeTracks) or as a group (e.g., with computeMatrix). If the peaks are obvious and you trust your peak calling then sure you can just use that instead, but given that that's prone to error it's usually better to simply use normalized signal.

Regarding actually looking at occupancy, I don't think that's useful. Occupancy isn't binary. If you're instead wanting to associate peaks with genes then using the nearest gene is common (it's not great, but it's not terrible either).
> 4 votes

---
Tags: visualization, chip-seq, peak-calling
---


title: Deep learning RNA sequences

text:
# Question
Title: Deep learning RNA sequences
Currently I'm working on a project, which combines deep learning with RNA sequences. I'll try to predict pseudotorsion angles \[1\] from raw rna sequence. The ideas is to train a neural network with raw rna sequences and for each nucleotide their corresponding pseudotorsion angles, and than to predict the angles of the remaining sequences in the test set.

How my data is structured: example:

```
seq1:   A  C  G    G   U   A   C
Eta:   169 87 110 87  45  187 78
Theta: 123 10 45  168 132 34  100

```

\[1\] These are angles describing the backbone conformation of a rna molecule.

I'm pretty new to the field of deep learning, and so far I build a simple feedforward neural network, but it's prediction accuracy is pretty low with only one percent.

Has anyone some tips for me how to improve this? How do I preprocess this kind of data correctly for deep learning?

I appreciate any help.

# Answer
> 2 votes
The question of folding of RNA sequences looks slightly similar to protein folding - perhaps searching in this domain might bring more suggestions.

An example of (current) state of the art of deep learning approach to protein folding is AlphaFold developed by DeepMind for CASP competition . They basically decompose the angles to an image representation and look for structures in the proximity of each atom. While it requires a **significant** processing power, RNA structures might be easier due to the size of the sequence.

# Answer
> 2 votes
I come from a protein background and this problem is analogous to protein torsion angle prediction, which in turn is a variant of protein secondary structure prediction.

Conventional ways to go here would be to use a sliding window of *x* bases either side to make a prediction for each point, or use a CNN/LSTM-type model to run over the whole sequence and output angles.

There are other technical concerns such as how to split the training and test sets, how to output the angles from the model (usually transformed by sin or cos), and whether to predict a single value or the probability of the value being in certain ranges.

The reference to protein tertiary prediction and CASP in another answer isn't directly relevant here unless you intend to build 3D models from your torsion angles, which you might very well do.

Some references on torsion angle prediction for proteins:

# Answer
> 0 votes
I think the solution is more in the DL.

I would look at simpler ML, using a random forest as a starting point and then maybe SVC. If you get some reasonable predictive power on this then move towards neural network/Tensorflow type calculation.

---
Tags: python, sequence-analysis, rna, machine-learning
---


# Question
Title: what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?
What is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets nowadays?

About 5 years or so ago, it used to be the latest release of the CSHL HapMap 12 populations that were part of the 1000 genomes project. These were:

CHB CHD MEX GIH TSI LWK ASW MKK HCB JPT CEU YRI

For example, the EnsEMBL project currently has these 12 populations available as a vcf here:

```
http://ftp.ensembl.org/pub/current_variation/vcf/homo_sapiens/

CSHL-HAPMAP-HAPMAP-CHB.vcf.gz
CSHL-HAPMAP-HAPMAP-CHD.vcf.gz
CSHL-HAPMAP-HAPMAP-MEX.vcf.gz
CSHL-HAPMAP-HAPMAP-GIH.vcf.gz
CSHL-HAPMAP-HAPMAP-TSI.vcf.gz
CSHL-HAPMAP-HAPMAP-LWK.vcf.gz
CSHL-HAPMAP-HAPMAP-ASW.vcf.gz
CSHL-HAPMAP-HAPMAP-MKK.vcf.gz
CSHL-HAPMAP-HapMap-HCB.vcf.gz
CSHL-HAPMAP-HapMap-JPT.vcf.gz
CSHL-HAPMAP-HapMap-CEU.vcf.gz
CSHL-HAPMAP-HapMap-YRI.vcf.gz

```

Each SNP has an AF entry from which a multi-populations vcf with rsids, alleles and frequencies can be built, where the AFs are such as AF\_CHB, AF\_CHD, etc.

The 1000 genomes project populations documentation describes more than these 12 populations, but I haven't seen equivalent population vcfs with AFs built from the individuals within the population for the remainder of these, apart from the original HapMap 12 marked with a `*` below:

```
###Populations and codes

*         CHB   Han Chinese             Han Chinese in Beijing, China
*         JPT   Japanese                Japanese in Tokyo, Japan
         CHS    Southern Han Chinese    Han Chinese South
         CDX    Dai Chinese             Chinese Dai in Xishuangbanna, China
         KHV    Kinh Vietnamese         Kinh in Ho Chi Minh City, Vietnam
*         CHD   Denver Chinese          Chinese in Denver, Colorado (pilot 3 only)

         CEU    CEPH                    Utah residents (CEPH) with Northern and Western European ancestry
*         TSI   Tuscan                  Toscani in Italia
         GBR    British                 British in England and Scotland
         FIN    Finnish                 Finnish in Finland
         IBS    Spanish                 Iberian populations in Spain

*         YRI   Yoruba                  Yoruba in Ibadan, Nigeria
*        LWK    Luhya                   Luhya in Webuye, Kenya
         GWD    Gambian                 Gambian in Western Division, The Gambia
         MSL    Mende                   Mende in Sierra Leone
         ESN    Esan                    Esan in Nigeria

*         ASW   African-American SW     African Ancestry in Southwest US
         ACB    African-Caribbean       African Caribbean in Barbados
         MXL    Mexican-American        Mexican Ancestry in Los Angeles, California
         PUR    Puerto Rican            Puerto Rican in Puerto Rico
         CLM    Colombian               Colombian in Medellin, Colombia
         PEL    Peruvian                Peruvian in Lima, Peru

*         GIH   Gujarati                Gujarati Indian in Houston, TX
         PJL    Punjabi                 Punjabi in Lahore, Pakistan
         BEB    Bengali                 Bengali in Bangladesh
         STU    Sri Lankan              Sri Lankan Tamil in the UK
         ITU    Indian                  Indian Telugu in the UK

```

I presume there is more public data nowadays that can be accessed to build a more complete vcf files of allele frequencies per population, with as many populations as possible to use as a reference dataset with new data.

Thanks in advance.

# Answer
> 4 votes
The best set I'm aware of can be downloaded from https://gnomad.broadinstitute.org .

---
Tags: vcf, genomics, 1000genomes, allele-frequency
---


title: Function to find ATP binding sites

text:
# Question
Title: Function to find ATP binding sites
I'm currently doing a python project where I need to write code that extracts the ATP binding site from a FASTA file using a regular expression.

I've tried the following code but it does not work. It opens a new file but the file is empty. I don't know what's wrong.

```
from Bio import SeqIO
import re
import os
path="C:/Users/nessr/Desktop/info projet"
os.chdir(path)
records=list(SeqIO.parse("plasmo.fasta","fasta"))
ATPM="[LIV]G[^P]G[^P][FYWMGSTNH][SGA][^PW][LIVCAT][^PD].[GSTACLIVMFY].{5,18}[LIVMFYWCSTAR][AIVP][LIVMFAGCKR]K"
patern=re.compile(r"(ATPM)")

def atp_binding(records,patern):

    for item in records:
        for m in patern.findall(str(item)):
            istart=m.start()
            iend=m.end()

            iend=iend+50

            print(istart,m.group())

file = open("output.txt","w")

for seq_record in SeqIO.parse(records, "fasta"):
    file.write(atp_binding(seq_record.records,patern)+"\n")
    file.write('-------------------\n')
file.close()

```

# Answer
The problem in your code is that the function atp\_binding is not returning any value, but printing it to the stdout.

You should add a return statement with the result you want to save in your output file.

Also, you should change how the resulting value is being saved. As the return value of your function is None, the file handler will write nothing.

I would use something similar to this:

```
def atp_binding(records,patern):
    s = ""
    for item in records:
        for m in patern.findall(str(item)):
            istart=m.start()
            iend=m.end()
            iend=iend+50
            s+="{0}\t{1}\n".format(istart,m.group())
    return s

file = open("output.txt","w")

for seq_record in SeqIO.parse(records, "fasta"):
   file.write(atp_binding(seq_record.records,patern)+"\n")
   file.write('-------------------\n')
file.close()

```
> 1 votes

---
Tags: fasta, phylogenetics, homework
---


title: How to use extractfeat?

text:
# Question
Title: How to use extractfeat?
I am working on a script on bash linux to get the CDS (coding sequence) of a gene using extractfeat by EMBOSS. It gives me the error: Warning: No sequences written to output file

I am unsure if my syntax is right, but this is what I have:

```
extractfeat inputFile outputFile -type CDS

```

# Answer
> 0 votes
Perhaps,removing introns (if it is a higher eukaryote)

```
extractfeat tembl:x65921 -type CDS -join

```

It is difficult to know without knowing the gene.

BTW NCBI's Entrez "efetch" will do essentially the same thing, but its more complicated

---
Tags: sequence-analysis, bash, emboss
---


# Question
Title: Sibling vs Parent/Child relationship detection using whole-genome sequence data
I am analyzing a family-based whole-genome sequence (WGS) dataset. My identity-by-descent(IBD) analysis shows they share about 50% of their DNA. I think my pedigree structure may have gotten mixed up though and I want to see if one individual is the mother or the sister of a male individual. Is there anyway to do this either with Plink or KING?

# Answer
With at least three samples, you can determine this by looking at the IBS0 output column (from either plink2 --make-king-table or KING): if the relationship in question has a very low IBS0 value compared to the other pairs, that's parent-child. You can try to do the same thing with only two samples, but with no other IBS0 values to compare against it's easier to misjudge.
> 1 votes

---
Tags: wgs, plink
---


# Question
Title: Is it possible to merge scRNAseq data from experiments with different design?
I have 4 different single-cell RNAseq experiments, each one representing a different sample of cell type population. I'd like to merge them to a single dataset. However, different cell types are enriched in each sample, and different samples can represent also different conditions (control and disease).

If I run a PCA on the merged dataset (without any processing of data) I find that the different sample's cells cluster separately. This could be due to biological and technical (i.e. batch effects) variation.

Is it possible to merge these different samples while retaining a good amount of information or should they be analyzed separately?

# Answer
It’s right there on the cellranger manual:

```
#aggregate results of counts for separate samples
cellranger aggr
#analyse the combined results
cellranger reanalyze

```

Note: I’m not sure this applies to the above case where batch effects could be an issue. However, this is how the analysis is performed for 10X data from multiple runs in principle.

cellranger aggr

cellranger reanalyze
> 2 votes

# Answer
See this paper from the Marioni group, where they propose a method for correcting batch effects between single cell sequencing experiments when each experiment contains different sub-populations:

> Large-scale single-cell RNA sequencing (scRNA-seq) data sets that are produced in different laboratories and at different times contain batch effects that may compromise the integration and interpretation of the data. Existing scRNA-seq analysis methods incorrectly assume that the composition of cell populations is either known or identical across batches. We present a strategy for batch correction based on the detection of mutual nearest neighbors (MNNs) in the high-dimensional expression space. Our approach does not rely on predefined or equal population compositions across batches; instead, it requires only that a subset of the population be shared between batches. We demonstrate the superiority of our approach compared with existing methods by using both simulated and real scRNA-seq data sets. Using multiple droplet-based scRNA-seq data sets, we demonstrate that our MNN batch-effect-correction method can be scaled to large numbers of cells
> 2 votes

# Answer
Have you tried seurat3? Here is the reference: https://satijalab.org/seurat/v3.0/pancreas\_integration\_label\_transfer.html
> 1 votes

---
Tags: rna-seq, scrnaseq, batch-effects, merge
---


# Question
Title: How to compare transcriptomic profiles of two cell types (single cell RNA-seq)?
I found this interesting Single RNA-seq data set in GEO, but I am not sure how to analyze it appropriately.

They have deposited transcriptomic profiles of human and mouse pancreatic islets (pancreatic cells: Beta cells, Delta, etc). The problem I see is that the different pancreatic cell types are not in equal numbers: 872 Beta cells and 214 Delta cells.

I want to see whether the expression of a set of genes are correlated between the two cell types. I don't think that I can visualize the correlation between the expression of the Beta and Delta cells on a scatter plot if they are unequal in number. Is there a workaround to include all the cells in the analysis, without losing good quality transcriptome of beta/delta cell?

What I would like to do is to compare the expression of two pancreatic cell types (Beta v.s. Delta cells) using a scatter plot. I have a set of genes that I am interested in, and would like to see their expression profile for beta and delta cells. As far as I know, it is only possible to visualize this using a scatterplot (X:Beta vs Y:Delta) where the number of isolated beta cells is the same as the number of delta cells.

**Question:** Given unequal number of isolated panacriatic cells, what would be an appropriate way to compare the expression profiles? Should I just ignore the extra ones?

Any idea how to approch this problem?

Thanks,

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

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

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

# Answer
Scatterplots are used to represent multiple dimensions of the same piece of data. A scatterplot would be used in this case to display *all* genes in your set on a single plot, but I would have some difficulty interpreting such a plot. If you want to look separately at each gene, then a density curve (or histogram) would be more appropriate.

Your problem is similar to common problems in flow cytometry, and I don't think there's any reason why similar approaches can't be used. Here's a histogram / density plot showing two different types of cells (negative control, and treatment condition), indicating the difference in fluorescence of a particular marker (somewhat analogous to gene expression) across all cells of each type. Note that the number of cells are different between the negative control and treatment:

This image was taken from here, which also demonstrates a 2D scatter plot (where each point is a *cell*, rather than a gene).
> 1 votes

# Answer
seurat3 may help. You can compare the marker genes from the 2 cell types or treat one cell type as a control. Please see the example analysis from: https://satijalab.org/seurat/v3.0/immune\_alignment.html
> 0 votes

---
Tags: python, transcriptome, scrnaseq
---


# Question
Title: How to extract values from second file on the basis common first column?
File1:

```
ENST00000000233 ARF5
ENST00000000412 M6PR
ENST00000001008 FKBP4
ENST00000001146 CYP26B1
ENST00000002501 DBNDD1
ENST00000002596 HS3ST1

```

File 2:

```
ENST00000000412 0.779142
ENST00000001008 0.738143
ENST00000002501 0.715315
ENST00000002596 0.713664

```

Desired Output

```
ENST00000000412 0.779142    M6PR
ENST00000001008 0.738143    FKBP4
ENST00000002501 0.715315    DBNDD1
ENST00000002596 0.713664    HS3ST1

```

# Answer
> 6 votes
I like to use the `join` command for this. For your example above, you can simply use the following:

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

```

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

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

```

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

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

```

# Answer
> 1 votes
Using awk, this can be done in the following way:

```
awk 'NR==FNR{a[$1]=$2; next}$1 in a{print $1, $2, a[$1]}' file1.txt file2.txt > same.txt

```

Explanation:

```
NR==FNR (along with the next command)

```

Will go through the first file only and create an associative array of of the geneIDs corresponding with the geneNames.

Then if the 1st column in the 2nd file is present in the array, a, it will print both columns from the file plus the geneName.

---
Tags: bash, awk
---


title: Comparing two genome annotations

text:
# Question
Title: Comparing two genome annotations
I have one-dimentional array (human genome). Also I have two annotations for it, we can think about them as different peaks (it's nucleosome and secondary structures). How can we find correlation and/or causation between this annotations? Currently we are forming hypothesis like "one element is always before other" and statistically testing them. Can we be smarter and use known tools or at least automate this?

# Answer
> 1 votes
You can express the annotations of the genome with a dot plot, x and y represent the two annotations along the genome. The you can summary the result with Spearman correlation coefficient(s) based on whole set of annotations or particular gene families.

# Answer
> 1 votes
I would start with metaplots for a visual/exploratory analysis to guide your statistical testing. Cut the genome into short segments e.g. centred on the centre of each nucleosome peak, overlay all the segments on top of each other, and plot the abundance of secondary structure features at each base count as you step away from the centres. This should give an easy view on whether there is any clear structure in the relative positions of the different objects.

Do the plots both ways round, and try different subsets of the objects (e.g. ones overlapping genes or not, different stringency thresholds for peak calling). This should hopefully give you some insights on what is worth testing for.

---
Tags: genome, statistics, sequence-annotation
---


# Question
Title: How to highlight specific cells in Seurat 2.4
I used Seurat 2.4 on our scRNA dataset to obtain the following tSNE plot.I was able to successfully extract cell IDs from the different clusters, and generate gene expression profiles. The analysis, and the biology makes sense. I have been however stuck in trying to highlight specific cells we are interested in using the Cell IDs (barcodes). For e.g, in this plot I am hoping to highlight specific cells (~ 120 of them) within Cluster 5 (blue cluster on the top). It can be either in featureplot mode or in this plot itself by an overlay, it doesn't matter. All I have to show are the 120 cells within the cluster. For eg. if cluster 5 has cells atag-1 atgc-2 atat-3, cacc-4 cat-5... i want to recreate this plot using atag-1 atgc-2, cacc-4

I tried using do.identify = TRUE but that is a visual selection and can select regions, but is not useful to select specific cells. I was thinking if there is a way to subset and replace cluster 5 to include just the cells (~120 odd out of the total 200) but was unable to crack it. I tried by using cells.use but then realized that doesn't work either.Any ideas will be greatly appreciated, to get me out of this rut. Thank you

# Answer
You can "highlight" any cells you'd like by creating a new column (or modifying an existing column) in the meta data slot within the Seurat object. From there you can modify the contents of that vector to your liking and plot the results.

It's not clear if you've been successful in obtaining the barcodes for the subpopulation to begin with, because that's a separate issue all together, and you need to clarify what attribute you are trying to subset the cells on because there are a number of ways to accomplish that. If you've already obtained the barcodes then I'll describe how you would plot the results below.

Let's assume you'd like to keep the original cluster IDs for all other cells. We'll start off by duplicating the cluster ID column (this isn't necessary, but it's good practice to preserve the original IDs).

```
n <- seurat_obj@meta.data # for convenience
n$new_column <- n$res.1.0

```

Now you need to modify the ids for cluster 5 to break them up into 2 populations.

```
n[rownames(n) %in% my_barcodes, "new_column"] <- "5b"

```

You can plot the new column using DimPlot with the argument `group.by = "new_column"`, which will make the subpopulation stand out from the rest of the cells in cluster 5.
> 1 votes

---
Tags: scrnaseq, seurat, visualization, single-cell, transcriptome
---


# Question
Title: Import gene list to Seurat to define cell types
Is there a way to import gene list into Seurat to define cell type? The default cell types in Seurat is not enough for our research. For example, we want to mark a subtype of B cells in Seurat, but seems like Seurat only have general B cell, I'm planning to download gene list from Genomic Cytometry, and then import the gene list into Seurat and then Seurat define cell types in the cluster based on the gene list I import. Does Seurat support this function?

# Answer
> 2 votes
You can use scID with your gene list to find the cells that match your gene list. You can then find which cluster in Seurat have these cells.

# Answer
> 1 votes
Seurat does not define cell types by name. It clusters and assigns each cell to a cluster, from 0 to X. If your data has the cell type (e.g. B,T, Mast cells) it means that someone annotate the clusters so that they have a biological meaning.

You can assign different names to the clusters by using the `AddMetaData` function.

```
sampleID  cellType
cell1     B
cell2     B
cell3     T
cell4     T

newMeta <- oldMEta
newMeta$origType <- newMeta$cellType
newMeta$cellType <- c("B1","B","Treg","Tkill")

object <- AddMetaData(object, newMeta)

```

The new metadata should look like

```
sampleID  origType  cellType
cell1     B         B1
cell2     B         B
cell3     T         Treg
cell4     T         Tkill

```

You can then swap the labels using `SetAllIdent(object = test, id = "cellType")`. Your tSNE or uMAp, will now show the new cell types.

Have a look at Manually define clusters in Seurat and determine marker genes

# Answer
> 0 votes
Maybe you can try `Seurat::AddModuleScore()`, then `FeaturePlot()` and see if some of your B cells are different. After plotting this on `GenePlot()`, perhaps you can set a cutoff, then assign identities.

Alternatively, use your B cell gene list in `RunPCA(object, pc.genes = yourgenelist)` instead of the usual variable genes. You can also subset B cells first with `SubsetData()`, and run the above.

---
Tags: r, seurat, single-cell
---


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

sample name:

```
188833984]$ ls

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

```

pair end file within sample name Leaf\_T1\_FD\_R10\_L001-ds.3b884c360b1e4ae185408a613b90a3bc

```
188833984]$ cd Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc

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

```

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

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

```
find . -name "*.fastq.gz"

```

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

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

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

```

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

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

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

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

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

```
> 3 votes

---
Tags: bash, awk
---


# Question
Title: extracting last few lines of a growing bgzip compressed .vcf.gz file
I am running `bcftools merge` to create a single `.vcf.gz` file out of a collection of chromosomal `.vcf.gz` files from the 1000g project.

The merge process has been running for a couple of days, and I want to estimate how much longer it will take. The file is a growing `.vcf.gz` that is being `bgzip -c` compressed from a Linux pipe:

```
bcftools merge [...] | some processing and filtering | bgzip -c > final.file.vcf.gz

```

How can I know at which point is the file currently, to predict how much longer it will take to complete the merge?

Is there a way to do something like:

```
tail --with-some-magic final.file.vcf.gz | bgzip -d -c

```

So that I can read at which chromosome is currently on?

Is there a way to guesstimate the final size of the merge file and compare from the sum of sizes of the original files? E.g. the sum of compressed sizes originally (before filtering and munging) is 20G:

```
find $PWD -name "ALL.chr*.vcf.gz" | xargs ls -l | awk '{print $5}' | csvtk summary -H -f 1:sum | numfmt --to-si

```

This gives 20G, but the merge file is at 27G and still running.

Any ideas?

# Answer
27GB is not huge. You can `zcat` the file and see the last few lines. It will take ~10 minutes. It is theoretically possible to peek a bgzf file, but no such tools are available.

By the way,

> I am running bcftools merge to create a single .vcf.gz file out of a collection of chromosomal .vcf.gz files from the 1000g project.

Probably you should use `bcftools concat`, not `bcftools merge`. `merge` combines VCFs from different samples. `concat` concatenates VCFs for individual chromosomes; each input VCF contains the same set of samples.
> 5 votes

# Answer
You could find which Chromosome the merged file is currently writing on by:

```
bcftools query -f'%CHROM\n' final.file.vcf.gz | tail

```

Note that also you might want to try multithreading the command that you are using to write the file, as bcftools is capable of using several threads for the compression of output files. E.g.:

```
bcftools merge -Ob --threads 32 file1.vcf.gz file2.vcf.gz > mergedFile.vcf.gz

```
> 0 votes

---
Tags: vcf, vcftools, bcftools
---


# Question
Title: creating population AF from set of individuals in 1000 genomes population?
What's the recommended way to create population-level Allele Frequency variant entries from a set of individuals (vcfs?) in the 1000 genomes (or similar) project?

For example, given the IBS individuals in the 1000 Genomes project:

http://www.internationalgenome.org/data-portal/population/IBS

I would like to create an IBS `.vcf.gz` that contains an IBS population-level AF entry for each of the variants present.

Any ideas on how to achieve this?

EDIT: this is related to my earlier question on how to create the most complete file for population allele frequencies, here: what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?

I have seen that there is an Allele Frequency Calculator tool in EnsEMBL, which is at least available with GRCh37:

http://www.internationalgenome.org/allele-frequency-calculator-documentation

http://grch37.ensembl.org/Homo\_sapiens/Tools/AlleleFrequency?db=core

There are also instructions on how to do this genome-wide using vcftools:

```
vcf-subset -c ~/1000g/IBS.samples.list ~/1000g/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | fill-an-ac | bgzip -c > ~/1000g/IBS.chr1.20130502.vcf.gz

```

This seems to be generating correct files, and the AF can be calculated with a tool that combines `AC` and `AN` values, such as it is described here:

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

Any other ideas / warnings / pitfalls of the approach?

# Answer
> 1 votes
bcftools has a plugin called fill-tags which should do the job for you easily, without having to mess around with lots of complicated bash scripting.

You could do something along the lines of this, presuming you have a .vcf file with all the IBS individuals in it:

```
bcftools +fill-tags in.vcf.gz -Oz -o out.vcf.gz

```

This will give each position in the .vcf file new INFO tags, one of which will be allele-frequency.

You can then pull off the frequency for a particular variant (CHR:POS) by:

```
bcftools query -t CHR:POS -f'%CHROM\t%POS\t%AF\n' out.vcf.gz

```

---
Tags: vcf, genomics, allele-frequency
---


# Question
Title: How to predict protein binding to ligands?
I want to predict from protein sequence if the protein binds to metal, nuclear or small ligand. How can I do this ? Which features are relevant if I want to use them in a machine learning algorithm ?

# Answer
> 1 votes
If you are training a machine learning algorithm you probably want to train against things you know to be true, rather than making predictions from the sequence and training against the predictions.

In this case you will want to use protein structures with the relevant ligands bound as your data. If you go to the RCSB PDB advanced search you can search for entries with any ligand present, specific ligands present, protein and DNA present, etc. Then you can extract the FASTA sequences for these entries.

See also databases like MetalPDB for specific ligand types, and one of our papers for an example of creating a dataset of metal binding proteins.

Features to use would be the raw sequence as a one-hot encoding on amino acid types. You could add other sequence level features such as secondary structure predictions, disorder predictions or physico-chemical properties. However a decent network would be able to predict much of this if it needs to. You could create alignments and put in a sequence profile or covariation data, which would be giving the network more information.

---
Tags: proteins, machine-learning, interactions, features
---


# Question
Title: Run Rscript for multiple samples (files)
I have used the following script for calculating the median for one file (gene). However, I have 20,000 files for which I want to run this R script and get 200 separate output files

```
 #!/usr/bin/R
 setwd("/location folder")

 d <- read.table("tgf.txt", header=FALSE)

 x <- aggregate(d[, 2], list(d$V3), median)

 write.table(x, "tgf.output.txt", quote=FALSE, row.names=TRUE)

 q()

```

# Answer
> 1 votes
Another option is to use `list.files` and `grep` to select only the `.txt` in your directory. A `for` loop can read-aggregate and save the files one at the time.

```
all.files <- list.files("/location folder")
txt.files <- grep(".txt",all.files,value=T)

for(i in txt.files){
  d <- read.table(paste("/location folder",i,sep=""),header=F)
  x <- aggregate(d[, 2], list(d$V3), median)
  write.table(x, paste0(i,".output.txt"), quote=FALSE, row.names=TRUE)
}

```

# Answer
> 3 votes
Please do not use `setwd` in a script, it changes working folder without explicit specifications from the user, so the results could be unpredictable for people that are not familiar with the content of the script.

Create a tab-separated file with input filename and output filename, or at least a file with all input files, unless you can pick them using a `list.files(pattern = "<pattern_here>")`. Read this file as a `data.frame` and use it with an `apply`, so a function processes the df row-by-row. Let's say the file (named `metadata.txt`) is:

```
in_file out_file
tgf.txt tgf.output.txt
abc.txt abc.filtered.txt
xyz.txt xyz.randomname.txt

```

Define a function called `process_row()` in the Rscript:

```
process_row <- function(r) {
    d <- read.table(r[[1]], header = FALSE);
    x <- aggregate(d[, 2], list(d$V3), median);
    write.table(x, r[[2]], quote=FALSE, row.names=TRUE);
}

```

And do the actual work:

```
meta_data <- read.table('./metadata.txt', sep = ' ', header = TRUE, stringsAsFactors = FALSE);
apply(meta_data, 1, process_row(r))
# You don't need a q() in an Rscript because it exits after the last line of code anyway.

```

---
Tags: r
---


# Question
Title: How to predict the distance between two residues in a protein sequence
I have some protein sequences and I need to predict the physical distance between each two residues in the protein when folded correctly in 3 dimensions; I need to predict this distances just by having the protein sequence and not the structure.

I have a group of pdb IDs and I have got homologs for them and made a multiple sequence alignment for each query sequence; But later I need to imagine that I don't have any structure for being able to predict the distance between residues in query sequences without any known structure.

There are some articles for solving this issue using machine learning methods but I don't want to train with a large numbers of structures. Any one knows any way to solve this problem by using mathematical or statistical methods?

# Answer
Developments in this area over the last year are worth pointing out. All use residue-residue covariation information from a protein family so require a multiple sequence alignment to work, though it sounds like you have that.

The most accurate way to do this would be to use pre-trained machine learning methods, which would not require you to train your own version. Examples include:

* DMPfold from our lab, which is freely-available open source.
* The RaptorX method from the Xu lab.
* The state of the art is AlphaFold from DeepMind, though this has not been published yet.

These predict the probabilities of a residue-residue distance being in various distance 'bins', e.g. the probability of being 5 - 6 Angstrom apart. They don't predict beyond ~20 Angstrom.

If you specifically wanted a statistical, non-trained method, then consider:

All of which are available in some form. The problem here is that they all predict binary contacts at 8 Angstrom, whereas it sounds like you need more fine-grained distances.

In a general sense, you can use any 3D structure prediction method and then read the distances back from the predicted structure.
> 1 votes

---
Tags: proteins, sequence-analysis
---


# Question
Title: swap genotypes for specific fields of a vcf file
Is there a tool already established to "swap" the ref/alt information in the `GT`, `AD` and `PL` fields of a `vcf` file?

For the `GT` field, I need to transform: `0/0 -> 1/1` ; `1/1 -> 0/0` ; `0/1 -> 1/0` ; and finally `1/0 -> 0/1`

For the `AD` field, I need to swap the columns: `x,y -> y,x`

For the `PL` field, I need to swap the columns: `x,y,z -> z,y,x`

Sample input:

```
Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  0/0:17,0:17:48:0,48,720 0/1:14,8:22:99:167,0,347 .....

```

Sample output:

```
    Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:720,48,0 1/0:8,14:22:99:347,0,167 .....

```

note: the same operations will need to be carried out for all samples if the VCF file has multiple samples (the `...` stand for other potential samples).

Before I dive into scripting ( basically playing with regex), I would like to know if something to process this is already available. I looked into bcftools and fixref but they don't seem to help for this.

# Answer
If we can always assume your fields will be `GT:AD:DP:GQ:PL`, you can use this relatively simple, brute approach awk script:

```
awk -vOFS="\t" -F"\t" '{
  ## just print the header lines
  if(/^#/){
    print;
  }
  else{
    ## The number of samples will be the extra number of fields
    ## in addition to the default VCF 9
    sampleNum = NF - 9
    ## Fix each sample
    for(i=1; i<=sampleNum; i++){
      ## get the fields
      num=split($(9+i), fields,":");
      ## genotype
      if(fields[1] ~ "0/0"){ fields[1] = "1/1"}
      else if(fields[1] == "1/1"){ fields[1] = "0/0" }
      else if(fields[1] == "0/1"){ fields[1] = "1/0" }
      else if(fields[1] == "1/0"){ fields[1] = "0/1" }
      ## AD
      split(fields[2],values,",")
      fields[2] = values[2]","values[1]
      ## PL
      split(fields[2],values,",")
      ## Now set these as the sample data
      $(9+i) = fields[1];
      for(k=2;k<=num;k++){
        $(9+i) = $(9+i)":"fields[k]
      }
    }
    print;
  }
}' file.vcf

```

On your example input, this returns:

```
Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:0,48,720 1/0:8,14:22:99:167,0,347    .....

```
> 1 votes

---
Tags: vcf, file-formats, variant-calling, text-processing
---


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

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

Operon structure of genome: \- https://github.com/philschatz/microbiology-book/blob/master/resources/OSC\_Microbio\_11\_07\_Operon.jpg

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# Answer
> 5 votes
> So, first of all, is my high level observation correct? Does the genome, at least for bacteria, have a global LISP like structure?

No. I'm afraid this is what it boils down to. I don't see any real similarity to LISP or other programming languages at all. Perhaps it would be clearer if you could explain *in detail* why you think they're similar, but as it stands, no, there is no such structure apart from a very simplistic similarity which is sort of the same as saying that the genome is like a car because it carries genes and could therefore be studied by a car mechanic.

# Answer
> 2 votes
You mention this as the base similarity you find:

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

As other has pointed, this just considers the genome as a sequence of where only the linear order is important. It doesn't take into account other things that happen on the genome that are know to affect the cells like the chromatin structure, the duplication of certain regions, the transversion, the effect of moving elements, methylation, or simply that a single gene can have multiple functions (pleiotropism, moonlighting) while a programming function does not.

Even if we considered it as code, certainly by using tools like HMM we could also learn more about the structure of the genes we have. While certainly, we need more tools than just the sequence to discover what does each gene do.

There are a lot more differences than similarities on this comparison, it can be helpful for certain aspects but it fails for other aspects of cells (without even going into the virus discussion or the mitochondrial).

# Answer
> 0 votes
The only possible parallel was the "processor 'core' wars" in Avida and Tierra. Tierra was initiated by Tom Ray, around 1990. Here they had a 'population' of compiled binaries were randomly mutated (yep the binary code was mutated) and reproduction was dependent on how much of the processor could be consumed. The simulation worked and an in silico ecosystem evolved, but was heavily dependent on the preconditions of the simulation and proved no more reliabe at understanding evolution than differential equestions. So it was a nice idea, and very cool it worked, but it came to zero.

That aside the best example was the concept of the "selfish gene", but has very little to do with programming languages IMO and everything to do with understanding kinselection.

# Answer
> -5 votes
Believe it or not, there are people who know more than you do about programming, and who also know a hell of a lot more than you do about the genome. This is not an avenue those people think is worth pursuing.

DNA is a chemical, and interacts with other chemicals based on its shape. There is no equivalent to that in a human designed programming language.

---
Tags: genome
---


# Question
Title: how to merge two fasta based on full fasta header and remove the duplicate
I am wondering if anyone can help to merge two fasta file based on the full header as initially gene ID may be similar in some of the sequences and to remove the duplicate entry based on full fasta header

I do have two multifasta file:

file1:

```
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

```

file2:

```
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

```

Expected output:

```
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

```

# Answer
You can use the `FastaToTbl` and `TblToFasta` scripts I have posted here for this:

```
FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa

```

If you run this on your example, you will get:

```
$ FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa
$ cat merged.fa
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATG
GTGCGGGAACGCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGT
TTGGATGTCTCTGGCGACGG
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGC
GAACGCTACGGGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGAT
GTCACCGACGACGGTGGGTG
>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAG
GTCGCGGTCGGCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTC
GACCGCCTCGAAGCCCGCCT
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAG
GTCGCGGTCGGCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTC
GACCGCCTCGAAGCCAGCGC
>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGC
CAGTTGGTACGGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGAC
GTCGGCCTCGATGTCACGGG

```

This will remove any cases where both the header and the sequence are identical. If you instead want to do this based on the header only, not taking the sequences into account, you can just process both files and only print sequences from the 2nd if their ID line hasn't been seen before. For example, with `awk`:

```
$ awk '{
        if(NR==FNR){
            if(/^>/){seen[$0]++;}
            print
        }
        else{
            if(/^>/){ skip=seen[$0] }
            if(!skip){ print }
        }
       }' file1.fa file2.fa
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

```
> 2 votes

---
Tags: fasta, text-processing
---


# Question
Title: Reducing a string to a certain amount of the first characters in R
So I have a list of several thousand RNA transcripts, and I am trying to find out how they are spliced. Right now I am running `matchLRpatterns()` from the Biostrings package with a max gap length of 0, after running a `matchpattern` function to categorize the transcripts by donor sites (where the first cut in an RNA transcript is made to cut out introns). But applying that to several thousand transcripts is quite time consuming, when you have 5 donor sequences to match to 9 acceptor sequences in 100 nucleotide long transcripts. Instead, I am thinking to sort by donor, remove all that comes before the splice junction by converting it to a string then using `gsub`, then using `matchpattern` again, but I only want to check the first 8 nucleotides after the splice junction for the consensus sequence that follows my list of acceptor sites. If I just run `matchpattern` without shortening the length of the transcript, I will just find acceptor sites in the RNA transcript, not the acceptor site at the splice junction.

I was thinking `lapply` the `gsub` function to the transcripts to remove up to the splice junction, then `lapply` the `strsplit` function to make a list of lists of groups of 1 character strings, `lapply` the `head` function to get the first 8 single character strings in each list, `lapply` the `paste` function to make them one string again, `lapply` the `DNAString` function to make the sequences readable by matchpattern, then lapply the `matchpattern` function to check the acceptor sequence at the splice junction. Is there a better way? This feels a bit clunky and roundabout way of doing it. I know gsub can remove characters after a certain character or sequence of characters, but I am not sure what the first sequence in each transcript is at the splice junction. Can I still use gsub, but tell it to remove characters after a certain number of characters? Or better yet, is there a way to get matchpattern to only check the first x nucleotides in a sequence?

# Answer
> 1 votes
Is there any particular reason why you want to do this via pattern matching in R? This sounds like something that SuperTranscripts and Lace might be useful for:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1284-1

> Finally, Lace will annotate each superTranscript with blocks either by using the graph structure itself or alternatively using a script called Mobius (packaged with Lace) which infers the annotation from spliced junctions discovered when mapping reads to the superTranscript.

---
Tags: r
---


# Question
Title: Python API for working with ENSEMBL genomes
Dear bioinformaticians,

I could find two python APIs to mine ENSEMBL genome databases. There is a `pyensembl` which is more popular and then there is `genomepy` which is less popular. Following is my experience with both.

**1. `pyensembl`:**

`pyensembl` is a fast and powerful API to mine ENSEMBL genome databases. I have been using it for a while and in my experience it works very well. The name of the package may make you believe that it is developed/maintained by ENSEMBL. That is not true. Still, being an open source package with a realistic potential to be the default API for mining ENSEMBL genome databases.

However, there are some issues related to it's development that users may find disconcerting.

For example, documentation. I could only find a comprehensive api documentation for release 0.8.10. Current release is 1.1. Not even release notes are available. Therefore users have to rely on few lines of documentation that is available in the README.md.

Also, as far as I know, it is not peer reviewed even though it is an excellent analysis resource.

**2. `genomepy`**

`genomepy` loads all of the genome files into the RAM, making it immediately unusable, especially if you are working with larger genome such as human genome.

So my question is: Is there any other better python API for working with ENSEMBL genomes?

# Answer
> 4 votes
The correct API for Ensembl is the Ensembl REST API which is updated and maintained by Ensembl, and language agnostic.

# Answer
> 6 votes
Now there's another option, the `ensembl_rest` module, a thin wrapper around the Ensembl REST API to simplify its usage and make it more pythonic. You can find the documentation here.

To clarify things, I'm the creator and maintainer, but still think it's a legitimate alternative.

---
Tags: python, genome, sequence-analysis, ensembl
---