# Question
Title: Identify non-coding regions from a genome annotation
I have this GTF file and I use the command below on a Linux machine to extract the coding regions of the genome:

```
awk '{if($3=="transcript" && $20=="\"protein_coding\";"){print $0}}' gencode.gtf

```

How I could do the inverse and keep only non coding regions?

I need non coding regions of human genome, but for asking my question here I referred to coding parts too.

After googling I see I actually need enhancers, promoters, or other regulatory elements of human genome to find cancer driver genes placed in these regions.

Likely non-coding regions of genome are here:

https://elifesciences.org/download/aHR0cHM6Ly9jZG4uZWxpZmVzY2llbmNlcy5vcmcvYXJ0aWNsZXMvMjE3NzgvZWxpZmUtMjE3Nzgtc3VwcDQtdjMueGxzeA==/elife-21778-supp4-v3.xlsx?\_hash=KQi5jfO3kT2c4Qw44j4Rg6YAyCBQilYuWHVYXcRDuuo%3D

# Answer
> 4 votes
Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR.

`UTR` has its own feature in the gtf file. So you can do this:

```
$ awk -v FS="\t" '$3=="UTR"' gencode.gtf

```

If the gtf file is compressed use this instead:

```
$ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"'

```

BTW: Why are you using such an old release of gencode? The current version is v29.

As @Wouter tells you, the non coding region of a genome is the complement of the coding regions. Coding regions have its own feature in the gtf file. You can get them with `$ awk -v FS="\t" '$3=="CDS"' gencode.gtf`. Reading the manual for bedtools complement is your task.

# Answer
> 3 votes
This isn't a problem that's easily solved with awk. It's not like you're extracting a feature that's annotated in the GTF file. Instead, you want the *empty space* between annotated features.

A few years ago I wrote a program called LocusPocus for a similar task. It uses a gene annotation to break down a genome into gene loci and intergenic regions. It handles overlapping annotations and other weirdness pretty robustly. The output will include both coding regions and non-coding regions, but you can identify the intergenic spaces as those with `iLocus_type` equal to `iiLocus` or `fiLocus`.

**Note**: the `--delta` parameter will extend each gene/transcript by 500bp by default.

**Caveat**: the program only accepts GFF3 input by default. Hopefully it won't be too hard to convert your GTF to GFF3.

**Another caveat**: eventual interpretation of these data will depend on what features are annotated in the genome and which annotations you include vs ignore. Do you want your non-coding regions to include non-coding genes, or should these be treated separately? Some non-coding regions will be full of transposable elements and other repetitive DNA, while others will have enhancers, promoters, or other regulatory elements. It's important to tread carefully before you jump to any conclusions.

# Answer
> 2 votes
If you want all transcripts from that gtf file whose type isn't "protein\_coding", you can use almost the same command, just change the `==` ("is") to `!=` ("isn't"):

```
awk '{if($3=="transcript" && $20!="\"protein_coding\";"){print $0}}' gencode.gtf

```

Or, a simpler version:

```
awk '$3=="transcript" && $20!="\"protein_coding\";"' gencode.gtf

```

Note that this will not include any of the havana transcripts in the file, but I am assuming that's what you want since that's what your original command did.

Specifically, the command will return the following types of transcript (the numbers on the left are the number of such transcripts in the file):

```
awk '$3=="transcript" && $20!="\"protein_coding\";"{print $20}' gencode.gtf  | sort | uniq -c | sort -nk1
      1 "translated_processed_pseudogene";
      2 "Mt_rRNA";
      3 "IG_J_pseudogene";
      3 "TR_D_gene";
      4 "TR_J_pseudogene";
      5 "TR_C_gene";
     10 "IG_C_pseudogene";
     18 "IG_C_gene";
     18 "IG_J_gene";
     22 "Mt_tRNA";
     25 "3prime_overlapping_ncrna";
     27 "TR_V_pseudogene";
     37 "IG_D_gene";
     58 "non_stop_decay";
     59 "polymorphic_pseudogene";
     74 "TR_J_gene";
     97 "TR_V_gene";
    144 "IG_V_gene";
    182 "unitary_pseudogene";
    196 "IG_V_pseudogene";
    330 "sense_overlapping";
    387 "pseudogene";
    442 "transcribed_processed_pseudogene";
    531 "rRNA";
    802 "sense_intronic";
    860 "transcribed_unprocessed_pseudogene";
   1529 "snoRNA";
   1923 "snRNA";
   2050 "misc_RNA";
   2549 "unprocessed_pseudogene";
   3116 "miRNA";
   9710 "antisense";
  10623 "processed_pseudogene";
  11780 "lincRNA";
  13052 "nonsense_mediated_decay";
  25955 "retained_intron";
  28082 "processed_transcript";

```

You might also want to remove that "translated\_processed\_pseudogene" since that is actually translated into protein and is therefore technically coding:

```
awk '$3=="transcript" &&
     $20!="\"protein_coding\";" &&
     $20!="\"translated_processed_pseudogene\";"' gencode.gtf

```

---
Tags: sequence-annotation, genome, gtf, text-processing, interval
---


# Question
Title: awk for motif splicing from an amino acid alignment
I'm looking to splice out specific amino acid alignment positions as a continuous motif from an alignment (~15 000 sequences). The code is written in Perl (I agree not many developers now) but splicing method is awk which a lot of you know. Please note the hash is "tied" (not shown), this just means the order of the sequences is preserved. I can get the awk code below to return single amino acid positions.

I am struggling to get a loop working in awk.

I agree the Perl wrapper makes it more complicated, but its part of a huge data pipeline (all written in Perl).

```
#!/usr/local/ActivePerl-5.20/bin/env perl
#use strict;
use warnings;

# I loaded the alignment into a hash, where key = sequence and value = protein
my %dna = (
    virus1 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF",
    virus2 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF",
    virus3 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF"
);

# I define the start and end of the alignment positions I am splicing
my $n = 1;
my $z = 5;

# I turn the alignment into a string with separated by carriage returns
my @array = map {"$dna{$_}"} keys %dna;
my $string = join "\n", @array;

# This is the awk loop, note in Perl the "``" means issue a shell command

my $splice = `echo "$string" | awk 'BEGIN {FS=""} {x = $n; while (x < $z); print \$ x; ++x}' 2>1`;
print $splice;
__DATA__

```

**Expected output**

```
LGVLV
LGVLV
LGVLV

```

**Observed output** I get the original sequences returned

I am familiar with Python/Perl solutions of line by line processing but awk is more efficient IMO.

# Answer
I *think* you are trying to use `awk` to get characters 3 through 7 (inclusive) of each input sequence. If so, a working version of your approach would be this:

```
my $splice = `echo "$string" | awk 'BEGIN {FS=\"\"} {x = $n; while (x < $z){ printf \"%s\",\$x; ++x} print ""}' 2>1`;

```

That, when added to your perl script, produces:

```
$ foo.pl
VLVIL
VLVIL
VLVIL

```

That isn't the output you show, but the output you show is just the first 5 letters of each sequence, which doesn't seem to correspond to what your script is doing.

Now, the reason your version didn't work, is that the `while` was like so:

```
while (x < $z); print \$ x; ++x}

```

Assuming the `\$ x` was a typo, the problem is that the `++x` isn't part of the while loop. So that should never even have exited. Seeing as we're in the same time zone, and given the time, I am guessing you're tired and that you miscopied something when trying to give us a reproducible example. In any case, to get the loop right, you need `{}`:

```
while (x < $z){print \$x; ++x}

```

However, if all you want is a substring, a loop is a very inefficient way of doing that. You can instead use `awk`'s `substr` command:

```
  substr(s, m[, n ])
             Return  the at most n-character substring of s that begins at position
             m, numbering from 1. If n is omitted, or if n specifies  more  charac‐
             ters than are left in the string, the length of the substring shall be
             limited by the length of the string s.

```

So, in your script, that would be:

```
my $splice = `echo "$string" | awk '{print substr(\$0,$n,5)}'

```

---

Since you're already using a language that excels in string manipulation, I really don't understand why you would try to use `awk` instead. It will not be faster or more efficient than the same thing in perl:

```
print substr($_,$n,5) . "\n" for @array;

```

By calling `awk` you are also calling the shell which will be slower by definition, especially if you're `echo`ing thousands of lines. In fact, I tried this with 15000 sequences (just repetitions of your first one) and using either the `awk` loop I gave above, or the `awk`+`substr`, or the native perl `substr`:

1. awk loop, average time taken to run the same thing 100 times: 0.0525
2. awk substr, average time taken to run the same thing 100 times: 0.0522s
3. perl substr, average time taken to run the same thing 100 times: 0.0464

And you can make it even faster by skipping the array and iterating over the keys directly:

```
print substr($dna{$_},$n,5) . "\n" for keys(%dna);

```

That averaged at 0.0376s over 100 runs.
> 1 votes

# Answer
If you have the sequences in a seperate file like this:

```
$ cat input
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF

```

you can use a simple `cut` to fetch a substring:

```
$ cut -c1-5 input
LGVLV
LGVLV
LGVLV

```
> 1 votes

---
Tags: perl, awk
---


# Question
Title: Why can I not install snakemake on my SLURM computer even though I can find it in the bioconda channel?
I searched for the package 'snakemake' on my SLURM cluster using:

`conda search --channel bioconda snakemake`

and I get many versions, up to 5.4.2.

I then try to install it using:

`conda install --channel bioconda --yes snakemake=5.4.2`

but it fails with the following message:

```
PackagesNotFoundError: The following packages are not available from current channels:

  - snakemake

Current channels:

  - https://repo.anaconda.com/pkgs/main/linux-64
  - https://repo.anaconda.com/pkgs/main/noarch
  - https://repo.anaconda.com/pkgs/free/linux-64
  - https://repo.anaconda.com/pkgs/free/noarch
  - https://repo.anaconda.com/pkgs/r/linux-64
  - https://repo.anaconda.com/pkgs/r/noarch
  - https://repo.anaconda.com/pkgs/pro/linux-64
  - https://repo.anaconda.com/pkgs/pro/noarch

```

Why?

Why is the bioconda channel not included in the current channels? I thought that was what the `--channel bioconda` did in my install command.

Note that `conda install --channel bioconda --yes samtools=1.9` worked fine and installed samtools successfully.

# Answer
I created a file .condarc in my home directory with the following content:

```
channels:
  - conda-forge
  - bioconda
  - defaults

```

and it works now. I don't understand why this is necessary though. I thought the channels in the .condarc file would be used if I did not specify them in my command.
> 2 votes

---
Tags: software-installation, snakemake, bioconda, anaconda
---


# Question
Title: Are there any Genome in a Bottle-like resources for non-humans (especially for invertebrates)?
Genome in a Bottle is an excellent resource that provides many types of DNA and RNA sequencing reads for a single individual/cell line to test genome assembly and analysis tools. For example there are Illumina WGS reads, Oxford Nanopore, PacBio, and 10X reads available for download at high coverage depth. In addition, there are pre-aligned bam files and verified variant calls.

Is there any aggregation of similar robust single-individual sequencing resources for non-human, non-mammalian metazoan species? The only way I can consider finding such resources at the moment is going through SRA projects one at a time and hoping for the best.

# Answer
In addition I would like to add "Vectorbase", the curator "Dan" is a good guy (I forget his surname right now) and it's a good team. However, please read my rather somber assessment of arthropod genomics first.

---

We appear to be talking about lower metazoans, e.g. not protists, or yeasts (the latter I suspect are really good for your question). I don't know anything about worms.

Its important to understand that arthropods are several (soccer) divisions lower than mammals in terms of the question you are asking, except for two model organisms.

Flybase, Drosophila melanogaster is the most important species for your purposes. The second most important species within the Class Insecta will be Anopheles gambiae (malaria mosquito).

Personally I would stick to the main players because overall whilst there are a lot of insect genomes (because there are a lot of insects), their genomes remain rough because of the difficulties in assembly and requirement to physical mapping their genomes in comparison to their genome size. Also except for D. melenogaster the (incredible) tardiness in producing a draft genome. For example the idea of PacBio is in its infancy for many entomologists: and this would help them a lot in assembling their genomes correctly. I agree things have picked a lot recently, but I remain a skeptic.

Within the arthropods I would be very surprised if the crustacea, myriapods or Chelicerata (e.g. spiders) were in better condition than hexapods.

If you are interested in protists certainly Plasmodium falciparum (African malaria) has the sort of information you seek (Plasmodb).

---

**Summary** I don't want to discourage you, insects are great, but you need to be selective if you are from a human/mammalian genetics background.
> 2 votes

---
Tags: phylogenetics, ngs, phylogeny, genomics, public-dataset
---


# Question
Title: programmatic secondary structure prediction for >36-mer DNA oligonucleotides
I'm writing a tool to automate Sanger sequencing primer design for a production lab that uses a universal-tail chemistry Sanger sequencing to verify NGS results. Essentially, the template DNA is amplification using primers with universal 5' tails (M13F/R), and then the "BigDye" reaction is performed with an M13F and M13R primers in separate wells.

I'm able to design logical primers for each region in a request using Primer3 (v2.5.0), but then I'd like to calculate folding thermodynamics for the M13-appended primers to determine the best primer partner for the 5' M13 appendages.

In other words, in the last step, I would like to know if I should use:

> M13F-forward\_primer + M13R-reverse\_primer

or

> M13R-forward\_primer + M13F-reverse\_primer

The problem is that Primer3 won't fold an oligo \> 36 nt, since the algorithm for Tm is diminishingly accurate for longer oligos. However, I believe thermodynamic predictions are valuable, even if only relatively, i.e. comparing one 5' appendage to another 5' appendage on the same allele-specific sequence, so I'd like the most painless answer, even if it isn't accurate in absolute terms.

I'm aware that UNAFold will do this, but it seems like extreme overkill for what I need to do.

For lack of a better solution, I've implemented a workaround wherein I pair each allele-specific primer with either M13F or M13R and submit the pair as a `check_primer` task. In other words, the left primer is one of the allele-specific primers that primer3 designed, and the right primer is either M13F, or M13R. I then choose the M13 appendage orientation based on the `PRIMER_PAIR_0_COMPL_ANY_TH` output values for the two possible orientations.

I would like to know: Is there a better solution? Primer3 docs say that the Tm calculation (and not necessarily folding thermodynamic calculations) is inaccurate for 37+ mers. Is there a way to force Primer3 to fold a sequence \> 36 nt? If not, is there a simple tool for folding \>36-mers?

**Follow-Up:** The workaround described in my OP was inadequate and UNAFold would cost $9k (which would take months for my employer to approve - don't get me started), so I tried @MichaelG's solution below. It turned out to work quite well. I'm computing the maximum Tm value for all fragments from every oligo, and that number agrees quite well with the browser-based UNAFold output.

# Answer
> 2 votes
If you really insisted on using Primer3 you could perform a 36bp sliding window analysis. You would submit your primers in 36bp 'chunks' with a step size of say 3 bp. You would need to "circularise" your primer within the sliding window to assess the "beginning" and the "end" and remove "artifactual" hairpins. Its doable via a script and you would automate the web upload and data extraction via a "digger/spider" or else just use a download equivalent (I've used a standalone). There would be a "blindspot" using this approach, which you could overcome if you wanted.

So summary, its complicated and I think your suggestion of using UNAFold is a better idea.

---
Tags: phylogenetics, phylogeny, software-recommendation, primer
---


# Question
Title: How can I remove (non-trivial) duplicates from a VCF file?
This is related to the question I asked here. Consider a vcf file that contains duplicate variants, but where the duplicates aren't simply the same thing in the same notation but instead one is a subset of the other. For example:

```
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

```

These the two variants are actually the same. They result in exactly the same genotype. Changing AACAC to AATAC at position 529514 just means change C to T at position 529516.

Is there any tool that can detect such duplicates and remove them? I tried `vcfuniq` from `vcflib`, but that doesn't seem to recognize this as a duplicate. I think it only looks at the 1st 4 fields and only considers duplicates those variants with exactly the same values in the 1st 4 fields:

```
$ ./bin/vcfuniq test.vcf
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

```

However, as explained in the linked question, EBI's `vcf_validator` considers this invalid. And it doesn't really make sense to have these duplicates in any case, so is there any way I can detect and remove them? Preferably an existing tool, but I am open to scripting solutions as well.

---

This is further complicated by cases like this one:

```
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  529514  AACAC   AAT,AATAC       0.00    .       .     GT    0/1
chr12   529516  529516  C       T       .       PASS    .     GT    0/1

```

Unfortunately, this one won't be caught by the approach in Daniel's clever script:

```
$ cat test2.vcf | foo.py
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  529514  AACAC   AAT,AATAC   0.00    .   .   GT  0/1
chr12   529516  529516  C   T   .   PASS    .   GT  0/1

```

# Answer
It turns out that `bcftools` can do this (tested on bcftools-1.8), if you give it the reference genome to test against:

```
$ bcftools norm -d none -f hg19.fa test.vcf
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -d none -f hg19.fa test.vcf; Date=Wed Feb 27 16:08:44 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  .       C       T       .       PASS    .       GT      0/1
Lines   total/split/realigned/skipped:  2/0/1/0

```

For the more complex case of the multi-allelic variant in the second VCF example from the question, you can run it through bcftools twice. Once using `norm` to left-align and split multi-allelic variants, and then again to remove the duplicates:

```
$ bcftools norm -m -any -NO z -O v -o - ~/test2.vcf |
    bcftools norm -d none -f hg19.fa
Lines   total/split/realigned/skipped:  2/1/0/0
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -m -any -NO z -O v -o - test2.vcf; Date=Wed Feb 27 18:18:32 2019
##bcftools_normCommand=norm -d none -f hg19.fa -; Date=Wed Feb 27 18:18:32 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  529514  CAC     T       0       .       .       GT      0/1
chr12   529516  529514  C       T       0       .       .       GT      0/0
Lines   total/split/realigned/skipped:  3/0/2/0

```
> 3 votes

# Answer
I'm no expert with VCF (few can say they are!) but I have worked a lot with VCF data in the last few years, both tools to consume and produce VCF. I've never seen variants encoded in this fashion, and it seems to be non-canonical. Typically:

* Single nucleotide variants (SNVs) are encoded with a single base as the REF allele and a single base as the ALT allele.
* In the case of insertions or deletions, the shorter of the REF and ALT alleles will be a single base, the base preceding the inserted/deleted sequence. Thus the first base of the REF and ALT alleles is always the same.
* In the rarer case of two or more consecutive substitutions forming a multinucleotide variant (MNV) the REF and ALT alleles will have the same length.

Using multi-bp strings of the same length to encode SNVs is unnecessary and, as you've pointed out, problematic. This makes me think its a bug or a "feature" of the variant predictor that produced the VCF.

In this case, I'd write a small script that would check for variants where the REF and ALT alleles have the same length. If the base is the same for REF and ALT in any position, drop it, and adjust the position accordingly.

The script below will convert these funky SNVs to the canonical representation, and will also work on MNVs. Standard tools should then work to remove the duplicates.

```
#!/usr/bin/env python3

def canonicalize(instream):
    for line in instream:
        if not line.startswith('#'):
            values = line.split('\t')
            pos = int(values[1])
            ref, alt = values[3:5]
            if len(ref) > 1 and len(ref) == len(alt):
                # How many bp to trim off the end
                for n, (r, a) in enumerate(zip(ref[::-1], alt[::-1])):
                    if r != a:
                        revoffset = -1 * n
                        break

                # How many bp to trim off the front
                for n, (r, a) in enumerate(zip(ref, alt)):
                    if r != a:
                        offset = n
                        values[1] = str(pos + offset)
                        values[3] = ref[offset:revoffset]
                        values[4] = alt[offset:revoffset]
                        break
                line = '\t'.join(values)
        yield line

if __name__ == '__main__':
    import sys
    for line in canonicalize(sys.stdin):
        print(line, end='')

```

**UPDATE**: Upon further reflection, the more complicated example you listed actually makes sense. In the reference we have the sequence `AACAC`, and the alternate alleles represent two variations on this: deletion of the final two bp (in the first case), and a point mutation of the middle `C` to `T` (in both cases). Usually, only a single bp precedes the definition of an indel, so I would've encoded this complex variant as `ref=ACAC alt=AT,ATAC`.

So the SNV is "implied by"/"encoded in"/"redundant with" the complex variant, but it's not strictly a duplicate. I'm curious whether the VCF validator complains about these cases as well?
> 2 votes

---
Tags: vcf, variation
---


# Question
Title: Are there any RepBase alternatives for genome-wide repeat element annotations?
I’m using the RepBase libraries in conjunction with RepeatMasker to get genome-wide repeat element annotations, in particular for transposable elements.

This works well enough, and seems to be the de facto standard in the field.

However, there are two issues with the use of RepBase, which is why I (and others) have been looking for alternatives (so far without success):

1. RepBase isn’t open data. Their academic license agreement includes a clause that *explicitly forbids dissemination of data derived from RepBase*. It’s unclear to what extent this is binding/enforceable, but it effectively prevents publishing at least some of the data I’m using and generating. This is unacceptable for open science.

* Subordinate to this, the subscription model of RepBase also makes it impossible to integrate RepBase into fully automated pipelines, because user interaction is required to subscribe to RepBase, and to provide the login credentials.
2. RepBase is heavily manually curated. This is both good and bad. Good, because manual curation of sequence data is often the most reliable form of curation. On the flip side, manual curation is inherently biased; and worse, it’s hard to quantify this bias — this is acknowledged by the RepBase maintainers.

# Answer
Dfam has recently launched a sister resource, **Dfam\_consensus**, whose stated aim is to replace RepBase. From the annoucement:

> Dfam\_consensus provides an open framework for the community to store both seed alignments (multiple alignments of instances for a given family) and the corresponding consensus sequence model.

Both RepeatMasker and RepeatModeler have been updated to support Dfam\_consensus.

I haven’t tried it yet but it looks promising.
> 14 votes

# Answer
For pre-existing reliabe TE libraries it is a bit of a mess, because not everybody deposits the species-specific TE libraries to a database like RepBase. And as far as I know DFAM contains only human resources, or am I wrong?

As for *de novo* generation of species-specific TE libraries (which should be done for any species not already present in eg. RepBase): There is no "gold-standard" how to tackle this best. In principle one has to think about two main parts -repeat detection -annotation

For **repeat detection** I would recommend using a combination of two things (which is necessary, because TE copies might miss in the assemblies as repetitive regions tend to be difficult to assemble and thrown away in the final assembly).

I) Repeat detection from raw reads (as with e.g. DNApipeTE or tedna or RepeatExplorer). For me, DNAPipeTE worked quite nicely, but everything has pros and cons. II) Repeat detection from assemblies (as with e.g. REPET or as mentioned before RepeatModeler)

Then the **annotation** of these repeats is tricky too, because most methods are relying on homology between the *de novo* TEs and the TEs from some (probably distantly) related species. But some programs take also structure into account (like REPCLASS). REPET can do both detection and annotation, but is a pain to get to run.

I would recommend using some programs to do de novo repeat detection on your species of interest on both the raw reads and assembly, clustering these libraries together (with e.g. uclust and 95% identity) and then run an annotation with homology and structural identification.

Probably the programs will **not give you complete, full-length TE**s but rather consensus sequences of several copies from TE families. If you want, you could search all copies of one family, extract them from the contigs plus boundaries and align them manually and curate boundaries manually. Then extend boundaries if not hitting the surrounding (non-alignable) regions or landmarks of TEs like LTRs or TIRs or so. But this is very time consuming if you only want to compare TE abundance between species for example, I would not do this and rather compare the abundance using read coverage (as in Bast et al. 2016). Depends all on the questions you want to ask.
> 7 votes

# Answer
You could use RepeatScout, which has defined repeat libraries for a limited number of species (including human, mouse, and rat). If your taxon is not represented, you can also do de novo repeat prediction with RepeatScout to build your own library to feed to RepeatMasker. The RepeatScout publication includes some comparisons with RepBase. Another related tool is RepeatModeler, which wraps RepeatScout with RECON and some other programs, and shares authors with the RepeatMasker team.

On the plus side RepeatScout/RepeatModeler are open source and do not use manual curation, meeting your criteria. On the negative, I'm not sure exactly how RepeatModeler and the component tools are maintained. The RepeatScout web and github pages have not been updated for several years, although the RepeatModeler page shows its latest release was in 2017. Anyway, I know that some combination of RepeatScout/RepeatModeler have been used to annotate repeats for some *fairly recent* newly sequenced genomes, e.g. for cichlids, coelacanth, and Darwin's finch, so I think it's fair to say this kind of approach is accepted in the field, at least for vertebrate genome projects.
> 5 votes

# Answer
AFAIK Dfam and Repbase are currently the two best sources of (a variety of) TE sequences.

In my genome annotations I have used RepeatModeler+RepeatMasker and then later used Repbase+tblastx and Dfam+nhmmer to classify them.

The classification process in my pipeline PhyLTR (https://github.com/mcsimenc/PhyLTR) is based on Dfam and Repbase. The process I used for LTR identification is

1. Putative ID with LTRHarvest (based on structural sequence characteristics)
2. Classification-by-homology to Repbase and Dfam
3. Removal of elements without homology to sequences in Repbase or Dfam.

This results in a set of LTR-Rs which are full-length and have evidence that they are LTR-Rs.
> 4 votes

# Answer
+1 for taking issue with RepBase.

I use the annotations from the Hammell Lab GTFs that they put out with TEtoolkit. It is similar to what you described to be using so this may be a redundant and useless answer but from the digging I’ve done they do seem to be comprehensive and well curated (for Drosophila, at least).
> 3 votes

# Answer
I know this question is a bit old, but this is still an issue for a lot of researchers not being able to access RepBase. It seems now that the most recent version of RepeatMasker is dependent upon RepBase for full functionality if masking anything other than human (currently DFAM only has human models). I recently discovered a *de novo* repeat masking approach called REpeat Detector (Red). This might be a solution for some looking to mask repeats on a genome assembly for annotation. The paper is here. I also then wrote a wrapper around Red to make it a bit more easy to soft-mask a genome, which you can find here.

One of the limitations with Red is that the repeats are not classified, so they are only identified. You would have to use some of the other tools mentioned above to try to classify them.
> 2 votes

---
Tags: sequence-annotation, database, repeat-elements, transposable-elements
---


# Question
Title: scRNASeq expression matrix with decimal values
I am trying to replicate some results of a scRNASeq experiment and, when I looked at the data provided by the author, I noticed that some of the counts in the expression matrix are represented as decimals.

Would this be because the authors uploaded a normalised count matrix or for other reasons? Also, please note that the authors used the 10X genomics `cellranger` pipeline for the processing of raw data.

I have provided a brief snippet of the data frame here:

```
structure(list(P1TLH_AAACCTGAGCAGCCTC_1 = c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), P1TLH_AAACCTGTCCTCATTA_1 = c(0,
0, 0.314759545319035, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.314759545319035, 0, 0), P1TLH_AAACCTGTCTAAGCCA_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1.18024827050334, 0, 0, 0, 0, 0, 1.18024827050334,
0, 0, 0), P1TLH_AAACGGGAGTAGGCCA_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.644905372424165, 0, 0, 0), P1TLH_AAACGGGGTTCGGGCT_1 = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.780522072779651,
0, 0, 0), P1TLH_AAAGCAACAGTAAGAT_1 = c(0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)), row.names = c("RP11-34P13.7",
"FO538757.2", "AP006222.2", "RP4-669L17.10", "RP5-857K21.4",
"RP11-206L10.9", "LINC00115", "FAM41C", "RP11-54O7.1", "RP11-54O7.3",
"SAMD11", "NOC2L", "KLHL17", "PLEKHN1", "RP11-54O7.17", "HES4",
"ISG15", "AGRN", "C1orf159", "RP11-465B22.8"), class = "data.frame")

```

Also, I am quite confused of the structure they used - as I don't have much experience with single cell datasets. Just to give some further context:

The `colnames` of the "full" dataframe are coded from P1TLH\_sequence\_1 to P5TLH\_sequence\_1. There are 5 patients in the study.

I'm assuming that each of the `colnames` represents a cell for a particular patient while the rows are the genes.

If this assumption is correct, how would I do the conversion from `colnames` in the dataframe to the "cell names"? And lastly, why are there decimals in this expression matrix?

Thank you all for your time,

# Answer
It looks like you have normalized counts. My guess is that they are log-normalized with some scale factor, probably 10000.

With the whole matrix you should be able to work out the scale factor for each cell. See https://www.nxn.se/valent/2018/10/25/unscaling-scaled-counts-in-scrna-seq-data

They typical format for scRNA-seq data is to have cells in the columns and genes in the rows. So here your rownames are the gene names and colnames are the cell names.
> 3 votes

---
Tags: r, scrnaseq, sequence-annotation
---


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

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

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

```

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

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

```

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

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

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

---
Tags: transcription-factors
---


# Question
Title: Converting VCF file to PLINK bed/bim/fam files
I am trying to find the best way to convert VCF files to PLINK binary bed/bim/fam files, but it seems like there are many varied ways to do this. (For example, using Plink 1.9 --vcf tag, bcftools, GATK, and vcftools).

Obviously they all probably have their uses and you could likely spend a while discussing this, but can anyone give a short, simple explanation as to whether the methods are generally the same, or whether one or another is generally best for certain situations? Thanks!

# Answer
Given that Plink reads in VCF files natively, and Plink is your desired output format, I expect that `plink --vcf <file>` would be the best option.

As mentioned by @sarah, Plink uses the sample IDs to create family/individual IDs, see "Variant Call Format" section on this page:

> VCF files just contain sample IDs, instead of the distinct family and within-family IDs tracked by PLINK. We offer three ways to convert these IDs:
>
> ```
> --double-id causes both family and within-family IDs to be set to the sample ID.
> --const-fid converts sample IDs to within-family IDs while setting all family IDs to a single value (default '0').
> --id-delim causes sample IDs to be parsed as <FID><delimiter><IID>; the default delimiter is '_'. If any sample ID
>
> ```
>
> does not contain exactly one instance of the delimiter, an error is normally reported; however, if you have simultaneously specified --double-id or --const-fid, PLINK will fall back on that approach to handle zero-delimiter IDs.
>
> If none of these three flags is present, the loader defaults to --double-id + --id-delim.

There are caveats associated with the types of variants that Plink will pull in from the VCF files, mostly related to limitations of Plink's own file formats. But if your end-goal is Plink, then that should be fine.
> 8 votes

# Answer
Additional tidbit I learned in the last year that I wanted to share with anyone else that's working on this.

Plink creates an extremely generic fam file for you, but if you are updating this fam file to reflect your population's family structure, sex, etc., note that when using plink binary files (bed/bim/fam) the order of the individuals in the fam file needs to match the order of the individuals in the other files (which is the order of the individuals in the generic fam file)

Do be careful if you care about which alleles are major/minor, as sometimes the order gets switched around in plink. For more info see "Variant Call Format" section on this page.
> 3 votes

---
Tags: vcf, plink, gatk, vcftools
---


# Question
Title: Plink 1.9 --merge-list changing order of fam file
I am doing some whole-genome sequence analysis and to speed things up split the data into chromosomes and also into regions within those chromosomes. (So I have my dataset split into ~300ish files) I am doing some IBD analysis with plink files and to do this IBD analysis I am combining the files into one big plink file.

I ran into issues today where when using plink 1.9 and the --bfile and --merge-list commands, the output fam file is different from the input fam file. My dataset has related individuals so the fam file is very important. Also its is important to note that when using binary files the order in the fam file needs to be exactly the same as the order the individuals appear in the other files (I learned this the hard way haha).

Does anyone know why plink 1.9 changes the order in the fam file when merging plink files, and how to fix this? It works with plink 1.07 but with 1.9 I get a different fam file.

Performing "diff" on input and output fam files shows they are different:

```
plink=$mydir/bin/plink_1_9/plink

$plink --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

```

Using the same code but using plink-1.07 creates matching fam files as expected"

```
plink1=bin/plink-1.07-x86_64/plink
$plink1 --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

```

The merge\_list.txt file is in this format:

```
$indir/chr22/286/$prefix.chr.22.region.286.qc.bed $indir/chr22/286/$prefix.chr.22.region.286.qc.bim $indir/chr22/286/$prefix.chr.22.region.286.qc.fam
$indir/chr22/287/$prefix.chr.22.region.287.qc.bed $indir/chr22/287/$prefix.chr.22.region.287.qc.bim $indir/chr22/287/$prefix.chr.22.region.287.qc.fam
$indir/chr22/288/$prefix.chr.22.region.288.qc.bed $indir/chr22/288/$prefix.chr.22.region.288.qc.bim $indir/chr22/288/$prefix.chr.22.region.288.qc.fam

```

Can someone explain to me why this could be happening and if there is a way to fix it? Thanks so much!!

(Note, I know I can manually replace the fam file with the correct one with cp or mv, but I would like to know how to handle this with plink 1.9, if possible, as I am afraid it could be indicative of an underlying error in the merge process :) )

# Answer
Thanks Christopher Chang for the great answer via the plink2 google group! (See here)

His answer was as follows:

"Yes, plink 1.9 did change the default merged sample order. However, you can request plink 1.07’s behavior by adding “--indiv-sort 0” during the merge. (You can also use --indiv-sort + --make-bed at any time to change the sample order to whatever you want.)

As for why the default was changed, take a look at the new default order: it will probably make a lot of sense. It’s “natural-sorted” which tends to make it easy for humans to manually search through regardless of what mix of letters and numbers are used in the IDs."

See also the Plink 1.9 documentation for --indiv-sort which makes the point that this behavior is default only for merges--in other commands Plink 1.9 uses the same ordering as in Plink 1.07, which is the order the samples were in.
> 2 votes

---
Tags: plink, wgs
---


# Question
Title: Why does my SplitDotPlotGG plot show differential expression of a gene but Seurat::FindMarkers() function does not return it?
I have analysed and clustered my single cell rna seq data with methods in the Seurat package. After clustering, I obtained two main clusters: 0 and 1. I used the Seurat **SplitDotPlotGG** function to visualise the differential expression of a gene of interest X.

```
SplitDotPlotGG(object,
               grouping.var = "clusterID",
               genes.plot = c("geneX"),
               cols.use = c("blue", "red"), plot.legend = T)

```

The plot shows clearly that there could be a significant differential expression of geneX between the conditions represented in clusters 0 and 1. However, I cannot see geneX in the list of differentially expressed genes returned when I run:

```
FindMarkers(object, ident.1 = 0, ident.2 = 1, genes.use = "geneX")

```

Is there be any explanation for this? Thanks in advance.

# Answer
> 2 votes
You're showing groups 1 and 2 in the plot but comparing groups 0 and 1 in the differential expression test, so you need to instead set `ident.1 = 1, ident.2 = 2` in `FindMarkers`.

The `FindMarkers` function by default only returns genes with \>0.25 log fold-change expression, with \>3 cells in each group. It's possible that these criteria are not being met for gene X. You can try changing these defaults, setting `logfc.threshold = 0` to return all genes.

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


# Question
Title: R studio failed to install the latest development version of hdf5r from Github (win10)
I install hdf5r following https://github.com/hhoeflin/hdf5r， it shows two methods: When I use the second method to install from github,it showed many errors：

```
The downloaded source packages are in
‘C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\downloaded_packages’
√ checking for file 'C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\remotes8c87faaf3d\hhoeflin-hdf5r-cc278c1/DESCRIPTION'

* preparing 'hdf5r': (4.9s)
  √  checking DESCRIPTION meta-information ...
* cleaning src
* checking for LF line-endings in source and make files and shell scripts (676ms)
* checking for empty or unneeded directories (1.5s)
* building 'hdf5r_1.0.1.tar.gz'
  Warning: file 'hdf5r/cleanup' did not have execute permissions: corrected
  Warning: file 'hdf5r/configure' did not have execute permissions: corrected

* installing _source_ package 'hdf5r' ...
  checking hdf5 headers and libraries
  Environment variable LIB_HDF5 is not set. Attempting to download HDF5 files from repository.
  试开URL’https://github.com/mannau/h5-libwin/archive/master.zip'
  Content type 'application/zip' length unknown
  ** libs
  错误: (由警告转换成)this package has a non-empty 'configure.win' file,
  so building only the main architecture
* removing 'D:/soft ware/R-3.5.2/library/hdf5r'
  In R CMD INSTALL
  Error in i.p(...) :
  (converted from warning) installation of package ‘C:/Users/dell/AppData/Local/Temp/RtmpI7TP9f/file8c863b6795/hdf5r_1.0.1.tar.gz’ had non-zero exit status
  In addition: Warning messages:
  1: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers
  2: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers

```

But when I use the first command:install.packages("hdf5r"), it worked:

Can anyone help me? I search many times, but I could't find the answer.

# Answer
> 2 votes
Try running this:

```
source("https://bioconductor.org/biocLite.R")
biocLite("hdf5r")

```

Just a different method of installing packages besides github. It is more centred around bioinformatics packages as well so it is pretty useful.

---
Tags: r, seurat, software-installation, rstudio
---


# Question
Title: In GFF3, annotating more than one protein-coding gene (i.e. polycistronic) contained in a eukaryotic mRNA
Below is an example of a simple GFF3 file:

```
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS 3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS 3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

```

My question is: if we found another coding sequence (encoding a different protein) range from 3752 to 3904, how should the GFF3 be formatted? It seems to me that the GFF3 file can only allow one protein-coding gene per mRNA. If not, could anyone show me one example? Thank you!

\[Edit\]: Sorry, using "CDS" is misleading. I changed the title and contents to "contain more than one protein-coding genes" instead of "more than one CDS". I also change the starting position of the second protein to 3752 (was 3750) to reflect the fact. My sincere apologies.

# Answer
> 5 votes
The GFF3 file format specification doesn't care where annotations are described in the file, although it can make some things easier if the annotations are ordered by the start position:

```
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3750    3904    .   -   0   ID=ATRV01011:CDS:1;Parent=ATRV01011.0
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

```

There are many examples in well-studied genomes of overlapping genes, so the GFF3 format needs to be flexible enough to handle these. Different genes can be encoded across the same region, but on different strands. But different genes can also be encoded on the same strand, sometimes sharing an exon.

# Answer
> 2 votes
It is certainly possible to encode polycistrons in GFF3. In fact, this is discussed explicitly in the GFF3 specification.

> ### *Polycistronic transcripts*
>
> This is the case in which a single (possibly spliced) transcript encodes multiple open reading frames that generate independent protein products.
>
> `----->XXXXXXX*-->BBBBBB*--->ZZZZ*-->AAAAAA*-----`
>
> Since the single transcript corresponds to multiple genes that can be identified by genetic analysis, the recommended solution here is to create four "gene" objects and make them the parent for a single transcript. The transcript will contain a single exon (in the unspliced case) and four separate CDSs:
>
> ```
> chrX  . gene XXXX YYYY  .  +  . ID=gene01;name=resA
> chrX  . gene XXXX YYYY  .  +  . ID=gene02;name=resB
> chrX  . gene XXXX YYYY  .  +  . ID=gene03;name=resX
> chrX  . gene XXXX YYYY  .  +  . ID=gene04;name=resZ
> chrX  . mRNA XXXX YYYY  .  +  . ID=tran01;Parent=gene01,gene02,gene03,gene04
> chrX  . exon XXXX YYYY  .  +  . ID=exon00001;Parent=tran01
> chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene01
> chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene02
> chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene03
> chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene04
>
> ```
>
> To disambiguate the relationship between which genes encode which CDSs, you may use the `Derives_from` relationship.

This example doesn't indicate what the coordinates of each feature should be, but I think it makes most sense to have the `gene` and `CDS` feature coordinates correspond to the region coding each protein product, and the `mRNA` and `exon` feature coordinates correspond to the entire span.

**BIG CAVEAT**: Polycistrons are not typically considered a normal feature of eukaryotic genomes—they're much more common in bacteria and archaea. However, there are some confirmed cases of eukaryotic polycistronic mRNAs. That said, if overlapping gene predictions are observed in a GFF3 file, one cannot immediately assume that these are polycistronic. One or both of the gene predictions might be technical artifacts, or real genes that are incorrectly annotated, or (more rarely) actual overlapping genes. A claim that multiple protein products are expressed from a single transcript would require some convincing experimental validation.

# Answer
> 1 votes
If you found another CDS from 3750 to 3904 then it would necessarily be on a different transcript, given that you already have a CDS from 3860 to 3913 and having overlapping CDSs on the same transcript doesn't make much biological sense (though I'm sure there's some weird organism out there that has a gene like this). So the likely GFF3 file would look like:

```
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.2;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010.2:exon:1;Parent=ATNG01010.2
1   T1  CDS     3759    3904    .   +   0   ID=ATNG01010.2:CDS:1;Parent=ATNG01010.2
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010.2:exon:2;Parent=ATNG01010.2
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010.2:CDS:2;Parent=ATNG01010.2
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010.2:exon:3;Parent=ATNG01010.2

```

Note that I've just added another transcript, since that's likely what would be biologically going on.

# Answer
> 1 votes
I found one similar example in an Arabidopsis annotation file. It is case of the cpuORF (CONSERVED PEPTIDE UPSTREAM OPEN READING FRAME ). Although the first and the third mRNA isoforms are identical, they repeated the mRNA isoform and name them differently to accommodate the additional ORF/CDS. I guess this is another way to do it.

See the ranges for AT5G03190.1 isoform and AT5G03190.3 isoform.

```
Chr5    TAIR10  gene    758374  760382  .   +   .   ID=AT5G03190;Note=protein_coding_gene;Name=AT5G03190
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.1;Parent=AT5G03190;Name=AT5G03190.1;Index=1
Chr5    TAIR10  protein 758793  760148  .   +   .   ID=AT5G03190.1-Protein;Name=AT5G03190.1;Derives_from=AT5G03190.1
Chr5    TAIR10  exon    758374  760382  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  five_prime_UTR  758374  758792  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  CDS 758793  760148  .   +   0   Parent=AT5G03190.1,AT5G03190.1-Protein;
Chr5    TAIR10  three_prime_UTR 760149  760382  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.2;Parent=AT5G03190;Name=AT5G03190.2;Index=1
Chr5    TAIR10  protein 758539  760148  .   +   .   ID=AT5G03190.2-Protein;Name=AT5G03190.2;Derives_from=AT5G03190.2
Chr5    TAIR10  exon    758374  758660  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  five_prime_UTR  758374  758538  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  CDS 758539  758660  .   +   0   Parent=AT5G03190.2,AT5G03190.2-Protein;
Chr5    TAIR10  exon    758843  760382  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  CDS 758843  760148  .   +   1   Parent=AT5G03190.2,AT5G03190.2-Protein;
Chr5    TAIR10  three_prime_UTR 760149  760382  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.3;Parent=AT5G03190;Name=AT5G03190.3;Index=1
Chr5    TAIR10  protein 758539  758676  .   +   .   ID=AT5G03190.3-Protein;Name=AT5G03190.3;Derives_from=AT5G03190.3
Chr5    TAIR10  exon    758374  760382  .   +   .   Parent=AT5G03190.3
Chr5    TAIR10  five_prime_UTR  758374  758538  .   +   .   Parent=AT5G03190.3
Chr5    TAIR10  CDS 758539  758676  .   +   0   Parent=AT5G03190.3,AT5G03190.3-Protein;
Chr5    TAIR10  three_prime_UTR 758677  760382  .   +   .   Parent=AT5G03190.3

```

---
Tags: sequence-annotation, gff3
---


# Question
Title: Normalization for two bulk RNA-Seq samples to enable reliable fold-change estimation between genes
I have two bulk RNA-Seq samples, already tpm-normalized.

I would like to know what is a reasonable normalization procedure to enable downstream log fold-change estimation.

The distribution of the two samples using the common set of genes looks similar:

However, the two samples have only been tpm-normalized, is it enough to guarantee reliable fold-change estimation? Should I use another normalization procedure, e.g. Quantile Normalization, before comparison?

My objective is to define a signature using the genes that are up-regulated in Sample1 with respect to Sample0, and vice versa. I'm using log fold-changes, but I'm concerned that their value may be affected by each sample distribution.

Do you also have suggestions for the definition of up-regulated genes with these data?

# Answer
> 5 votes
You have only two samples?

You aren't going to be able to draw strong conclusions from that no matter what you do. Clever statistics don't work without replicates.

# Answer
> 3 votes
It's not a good idea to do tpm normalisation prior to differential expression analysis, because the actual read counts are useful to determine shot noise and statistical significance. DESeq2 includes read normalisation as part of its methods for differential expression analysis.

I think shot noise is best explained in terms of shooting photons at a target, as found on Wikipedia:

> Shot noise exists because phenomena such as light and electric current consist of the movement of discrete (also called "quantized") 'packets'. Consider light—a stream of discrete photons—coming out of a laser pointer and hitting a wall to create a visible spot. The fundamental physical processes that govern light emission are such that these photons are emitted from the laser at random times; but the many billions of photons needed to create a spot are so many that the brightness, the number of photons per unit of time, varies only infinitesimally with time. However, if the laser brightness is reduced until only a handful of photons hit the wall every second, the relative fluctuations in number of photons, i.e., brightness, will be significant, just as when tossing a coin a few times. These fluctuations are shot noise.

Sequencing also has shot noise effects because the sequencing process is a random sampling method. For low-abundance targets, the likelihood of sampling a target is low, so the effect of a single hit is amplified. This reduces the statistical significance when comparing different low-abundance expression values, because the variation in expression is highly dependent on which targets won the sampling lottery. The shot noise variation / expression curve can be seen in Figure 7 of this paper.

# Answer
> 3 votes
What I have generally done in the past is to process the data using `voom` in the `limma` package for bulk RNASeq. Inside `voom` you can call for different normalization methods to be used - "TMM" works fine for me and, is advocated by many in the field.

`voom` will output an object containing the normalized expression values in a log2 scale, which, also in my experience, has worked out just fine for calculating log fold changes.

Check out this link for more info on the package as well as normalization methods: https://www.bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html It is a very thorough introduction to the package and all of its capabilities.

Good luck!

---
Tags: rna-seq, normalization, fold-change
---


# Question
Title: Why are the clusters in the left side of a tSNE plot and not in the centre?
I have no idea why I used the subsetdata in Seurat after using CCA the TSNEPlot shows the cluster all in the left side and merge， did I miss something? Here are the code and figure I got.

```
cluster <-SubsetData(object = immune.combined, ident.use= "3")

cluster <- NormalizeData(object = cluster, normalization.method = "logNormalize", scale.factor = 10000)

cluster <- FindVariableGenes(object = cluster, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0, x.high.cutoff = 6, y.cutoff = 0, do.plot = F)

cluster <- ScaleData(object = cluster, display.progress = FALSE)
cluster <- RunPCA(object = cluster, pc.genes = cluster@var.genes, pcs.print = 1:10, genes.print =5)

PCAPlot(object = cluster, group.by= "orig.ident", dim.1 = 1, dim.2 = 2)

cluster <- FindClusters (object = cluster, reduction.type = "pca", dims.use = 1:10, resolution =0.3, print.output = 0, save.SNN = TRUE, force.recalc = T)

TSNEPlot(object = cluster, do.label = TRUE,do.return = TRUE,pt.size = 2,label.size = 5, no.legend=FALSE)

```

# Answer
> 1 votes
The issue here is simply that you did not re-compute the tSNE after subsetting the data, so you're looking at a subset of cells from a larger embedding. You would never get an embedding that looks like that by running tSNE on only those cells.

You also don't need to re-normalize the data after subsetting the cells. Since this is a per-cell calculation, the results will be identical to those in the larger matrix (you can check this for yourself).

# Answer
> 1 votes
I'm going to stick my neck out here. Overall PCA followed by tSNE is a very powerful method of unsupervised learning however, I agree with @GWW,

```
cluster <- ScaleData(object = cluster, display.progress = FALSE)

```

If this is a max-min rescale transformation (common in supervised learning) then the outliers will badly affect the clustering. You can remove outliers that are over 2 sd from a formal analysis and re-run the calculation. That is just the risk of a max-min rescaling.

The one issue I'm concerned about is the initial quasi-standardisation followed by log-transformation,

```
log1p(value/colSums[cell-idx] *scale_factor)

```

In this case you are multiplying by 10000, this is to remove the ratio (used as a quasi-standardisation) and then taking the log ??? That is not the way to transform a ratio, it would be better to take the reciprocal of the ratio. There are other (better) methods, a formal standardisation would be useful, before a transformation is considered. Possibly if that was done correctly you wouldn't have the problem of outliers potentially crunching your data.

If you fiddle with the parameters you might get something - but its a symptom not the cause of the problem. However, if removing the outliers works (I have not checked the ScaleData function), then you are sorted and all else is forgotten (at least for this data set). The information I read was here.

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


title: Understanding PCHeatmap outputs

text:
# Question
Title: Understanding PCHeatmap outputs
I am currently trying to understand the purpose of these PCHeatmaps - part of the `seurat` package in R:

All the online documentation I have searched for has only highlighted how it is used to observe heterogeneity in the data.. This much I understand.

My questions are:

* Are the genes in the rows of each PC plot just the most variable x number of genes defined in `FindVariableGenes`?
* If not, then how are each of these genes 'chosen' for use in the `PCHeatmap` visualization?
* What would be the different conclusions/deductions drawn from the PC1 heatmap and PC2 heatmap with regards to heterogeneity for example?
* And finally, how can these conclusions help in the analysis of scRNASeq data?

# Answer
> 3 votes
The `PCHeatmap` function (renamed `DimHeatmap` in Seurat v3) can be used to help determine the number of principal components to use in downstream analysis, as well as to visualize the top genes contributing to each PC.

Both cells and genes are ordered by their PC scores, and by default the 15 genes with the highest and 15 genes with the lowest PC loadings are displayed for each PC.

In your case with PC1 and PC2, my conclusion would be that PC1 and PC2 both strongly separate different populations of cells. For PC1 there is a smaller population with very strong separation, and PC2 splits the cells roughly in half.

You can also set the `cells.use` parameter (renamed `cells` in Seurat v3) to only show the top *n* cells for each PC which can help with visualization.

---
Tags: r, scrnaseq, statistics, seurat, pca
---


# Question
Title: Is nanopolish worth it since faster polishing software is available?
For Oxford Nanopore contigs produced by any long-read assembler has anyone performed any benchmarks to compare the polishing tools `racon`, `nanopolish`, or something like `wtpoa-cns`?

The compute time of `nanopolish` is quite long and I would rather use another tool if it will suffice. After long-read polishing I will polish more with two rounds of `pilon` using 200x coverage Illumina shotgun data from the same individual.

# Answer
> 3 votes
Nanopolish is necessary if you want to get high-quality consensus. racon etc don't use signal data. They can't achieve high quality. pilon at times doesn't work well. Even if it worked perfectly, it wouldn't help repetitive regions inaccessible by short reads.

---
Tags: nanopore, assembly, long-reads
---


# Question
Title: Repetition-free Longest Common Subsequence as similarity measure
I'm starting a thesis in Computer Science with topic the RFLCS (Repetition-free Longest Common Subsequence). I know very little about Bioinformatics problems and Genomics. I would like to understand the motivations behind this problem. More precisely, I'd like to understand why one need a similarity measure between genomes, and how the RFLCS is an useful measure.

\[Adi et al.\] introduce the problem in this article.

# Answer
A brief literature review can give you the answers.

This preprint gives the definition of RFLCS more clearly:

> A repetition free Longest Common Subsequence (LCS) of two sequences x and y is an LCS of x and y where each symbol may appear at most once.

Adi et al proposed to solve the genome rearrangement problem raised in this paper. If you want to know more, follow the references and citations of these papers and read them.

Basically, Adi et al think the length of RFLCS measures the distance between two genomes. The longer the more similar. I am not sure if this is a good measurement, though.
> 1 votes

---
Tags: genomics
---


# Question
Title: scater: SingleCellExperiment, Error in seq_len(ncol(assay))
I am trying to go through the following scRNA-seq tutorial.

But the line `sce <- newSCESet(countData=all.counts)` is not working anymore with the most up-to-date version of `scater`. Now we should use `SingleCellExperiment` function instead, so when trying to do the following:

```
sce <- SingleCellExperiment(assays = list(counts = all.counts))

```

I am getting the error:

```
Error in seq_len(ncol(assay)) :
  argument must be coercible to non-negative integer
In addition: Warning message:
In seq_len(ncol(assay)) : first element used of 'length.out' argument

```

Any help would be greatly appreciated.

# Answer
> 5 votes
I think you might want to use the most recent workflow, which should be for the latest version on Bioconductor.

# Answer
> 1 votes
You simply need to transform your data.frame into a matrix. Just use these two lines:

```
library(SingleCellExperiment)
sce <- SingleCellExperiment(assays = list(counts = as.matrix(all.counts)))

```

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


title: Significance and timing of "mux scans"

text:
# Question
Title: Significance and timing of "mux scans"
I'm using MinIONQC to do quality control on some ONT data. The software plots several read characteristics over time (hours passed during the sequencing process). These plots contain several vertical red lines. From the documentation:

> Muxes, which occur every 8 hours, are shown as red dashed lines

Reading this paper (doi:10.12688/f1000research.7201.1) I was able to get a rough idea of what a "mux" or a "mux scan" is. If I understand correctly, it's a calibration step that scans each well containing 4 nanopores to determine which nanopores will be used for data collection in which order.

Do I understand the basic idea of the mux scan correctly? Does it occur at the same frequency on all ONT instruments? Does it always occur in 8hr intervals, even if the sequencing reaction is run for longer than the standard run time?

# Answer
Yes, your understanding is largely correct. This originates from the situation that for each detector on a nanopore array there are 4 pores. I'll explain mux scans and groups, but this is outdated information as now another system is used.

So a MinION FC has 2048 pores with 512 sensors. Not all of these pores will be equally suitable for sequencing. At the beginning of the run originally a mux scan was performed, to test for each sensor which the best pore would be (if any) out of the 4 available. The best pores will be used for sequencing first, the rest will be turned off. That is group 1. After 8 hours the expectation is that many of these pores have degraded/accumulated damage or are not functional for other reasons. Then it is time to swith to group 2. So you change to another set of pores, which in these cases often lead to a temporary increase in yield. The default script used 8 hours, although some used adaptive software to select new pores when the first group started going bad. Also note that during the run no new mux scan was performed: the software did not check if the selection it made in the beginning still made sense. Pores in group 1 got turned off after 8 hours and switched to group 2, without taking their current 'health' into account. Cumulative yield plots get a characteristic bumpy pattern from this. Each bump is less than the previous one, eventually flatting off.

Now that's history, as the software now uses dynamic muxing and will continuously check for the best pore per channel - and as such maximize yield by having as many healthy pores as possible available for sequencing. Cumulative yield plots are now smoother.
> 6 votes

---
Tags: ngs, nanopore, quality-control
---


# Question
Title: Why does a missing label 11plex TMT shows up at almost 50% intensity compared to other labels?
We have a test plate that was done with 10 of the 11plex TMT labels. So, the intensities of the 11th label should be almost zero, plus the leakage from the 10th label due to 13C content. Although we observe a substantial drop in the intensities, we see ratios up to 50 percent compared to the highest reference intensity, which is the maximum of the intensities of the other labels, if we sum up the intensities for all proteins. For individual proteins the intensity sometimes is higher than for all "real" samples. Only talking about proteins that are identified "Confidential". Is that something you would expect? Is this due to the normalization?

I was using Reporter 0.7.17 from command line without specifying a reference sample. Pretty much default configuration and based on the default protein report.

Fig: Ratios (intensities) reported by the program.

# Answer
The developers answered along these lines: This might be a normalization artefact. Reporter is normalizing the ratios so that the median is 1. So, if the label is almost not present the noise gets increased. To check that you can turn of the normalization. This function is currently not documented. The developers are working on it though.

GitHub issue
> 1 votes

---
Tags: proteomics, tmt
---


# Question
Title: how to write the script to know each sample's cluster number?
I have followed the Seurat Tutorial:Integrating stimulated vs. control PBMC datasets to learn cell-type specific responses by using my data. And I am confused about how to write the script to calculate the cell number of each sample on each cluster.

I have run TSENPlot to know the samples and I'm using another TSENPlot to get the clusters. If I want to know find the cluster1 of each sample's number how can I write the script?

# Answer
Assuming you are following the basic seurat guide:

```
pbmc@ident

```

You can access the `'ident'` slot of your Seurat object and you will get a vector of the cluster IDs and the names will be the respective cells.

```
head(seuratobj@ident)
Cell1A7 Cell1C4 Cell1D1 Cell1D2 Cell1D3 Cell1D5
0 3 3 1 3 1
Levels: 0 1 2 3 4

```
> 2 votes

---
Tags: single-cell
---


title: Empty .result in MirDeep Star

text:
# Question
Title: Empty .result in MirDeep Star
I am using MirDeep star tool for RNA sequence alignment. I have tried both GUI version and command line version of this tool. But every time, I am getting empty . result file and .known\_mir file. This tool is neither giving any error nor any warning with my data. I have also tried the demo file given with this and this tool is working fine with this (i.e. giving some results in .result and known\_mir file).

Can anyone suggest me what could be the possible reason for this. Why I am getting empty file.

# Answer
I have found reason for not getting any result. I was using default adapter given by mirdeep\* which gives me zero alignment with my data. This error pointed out by author of mirdeep\* ( thanks to Jiyuan and colleen). After giving right adapter, I got the right result.
> 0 votes

---
Tags: rna-seq, sequence-alignment
---


# Question
Title: Unable to install bedtools on windows 10 ubuntu
I am trying to install bedtools on windows 10, but I get an error I don't understand. How can I fix it?

```
Building BEDTools:
=========================================================
DETECTED_VERSION = v2.27.1
CURRENT_VERSION  =
Updating version file.
 * Creating BamTools API
- Building in src/utils/FileRecordTools
  * compiling FileRecordMgr.cpp
make[1]: g++: Command not found
Makefile:36: recipe for target '../../../obj//FileRecordMgr.o' failed
make[1]: *** [../../../obj//FileRecordMgr.o] Error 127
Makefile:186: recipe for target 'src/utils/FileRecordTools' failed
make: *** [src/utils/FileRecordTools] Error 2

```

# Answer
> 3 votes
My personal advice would be to avoid compiling tools like this if you can avoid it, especially where you are using a "barebones" setup like the one you get with WSL Ubunutu as you'll almost always run into dependency issues.

The linux version of `conda` installs just fine in WSL and `bioconda` hae pretty much any bioinformatics package you might need, already compiled.

# Answer
> 2 votes
> make: g++: Command not found

You do not have g++ installed on your machine. This is used to compile the bedtools into an executable file.

Here is a tutorial on installing g++

---
Tags: genome, software-installation, bedtools
---


# Question
Title: Where can I find .cdf files for Affymetrix ChIP-chip data analysis?
I am analysing microarray ChIP data generated several years ago. I want to use the *affy* package for R. However, in order to load the CEL files and proceed with the analysis, I need the pombe annotation file which comes under the name of **Sp20b\_M\_v04.cdf**. I have looked everywhere for this file but I am unable to locate it. Is anyone in possession of this file or know where I can download it from?

# Answer
It looks like it's in GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL10187

There is apparently a zipped copy linked at the bottom in the supplementary files.

When people upload a study based on an array, there is a corresponding 'platform' record where the array details should be deposited.

Edit: you will probably need to gunzip the file first, to get it into CDF format.
> 1 votes

---
Tags: r, microarray
---


title: GTF from consensus sequence

text:
# Question
Title: GTF from consensus sequence
I am new to using bioinformatic tools and I was hoping this community could help clear something up for me. I need to generate a gtf file. My data are a set of complete HA genes for influenza B viruses in fasta format. From reading through forums and other internet searches, my understanding is that fasta cannot be transformed to gtf or gff because these last two provide location information which fasta does not (unless available in the header). I saw someone suggest mapping sequences to a reference sequence and generating the gtf from that mapping.

My question is, would it make sense to use the consensus sequence selected for alignment as this reference sequence? The alignment algorithm is a simple majority rule. Because all my data are complete genes, I would think the mapping of locations (start codon and stop codon) would be the same no matter which sequence in my data set I choose to use. Is this not so?

I am looking at viral evolution in the HA segment and my end goal is to calculate the dS/dN ratio for a set influenza B viruses collected during a particular season. The program I want to use requires a gtf, but likely other options would only require the reference genome. Would it make more sense to use the strain used as the vaccine candidate as the reference genome?

I appreciate any help you can provide me!

# Answer
> 1 votes
I was wanting more time, but okay here's what I'm thinking,

> Would it make more sense to use the strain used as the vaccine candidate as the reference genome?

Neither, because the GTF format is based on an HA protein structure (likely a crystal structure but could be cryoEM). This is not the same as an alignment position within a surface antigen, because most alignments of pathogen surface antigen genes will contain indels, e.g. between the GTF and your influenza epidemic. If HA doesn't have any indels fine, but I would be surprised.

The indels will cause artifactual shifts in the GTF domain boundaries unless they are corrected. Put simply an alignment position (particularly for pathogen surface antigen proteins) is rarely the same as the amino acid position of GFT because a protein structure of a surface antigen doesn't comprise indels in its native form.

If you do this work in detail, i.e. shifting between protein structure and phylogeny, you always write scripts which convert between aligned positions and "reference" (crystal/ cryoEM structure) positions. Its a common way of thinking.

So in summary, your reference sequence is the GTF sequence because that is the basis of the structure\<-\> alignment conversion

---

For positive selection analysis it does not matter what the reference is, **all that matters is the correct rooting**. For demarcation of domains within HA the GFT sequence needs to be at the top of the alignment 1) because it helps easily identify whether it contains indels 2) structure programs will expect the reference sequence to be there. If you are using the vaccine in a positive selection analysis ... you need to think hard about the question you are asking and the coverage of the vaccine. Epidemiologically speaking thats complicated because the passage of the vaccine will be an important source of selection and you will mistake that as the wild-type "escaping" the vaccine.

# Answer
> 3 votes
"Transforming" Fasta files to GTF or GFF3 files is a common request, as are similar tasks such as "converting" BAM files to VCF files. But these are underdefined questions, and "transform" or "convert" isn't really the right way to think about the task.

As far as your particular request is concerned, Fasta files and GTF files contain fundamentally different types data.

* Fasta files contain sequences, which could be chromosomes, assembled contigs, polypeptides, transcripts, or other biological sequences.
* GTF (and GFF3 and BED) files contain *intervals*, indicating that some feature of interest occurs at such-and-such position on such-and-such sequence. A hypothetical example: a transcription factor binding site occurs on chromosome 2 from position 40011 to 40017.

So you are correct: it's not really helpful to think of your task in terms of "transforming" Fasta to GTF. The process of creating a GTF/GFF3/BED file is often called *annotation*. There are many different methods to annotate many different types of genomic features: transposable elements, segmental duplications, protein-coding genes, non-coding RNAs, and so on.

As a result, it's meaningless to say "I need a GTF file" without clearly specifying what that GTF file will annotate. (It's kind of like saying "I need an XML file" or "I need an Excel spreadsheet" without indicating what the contents of those files would be.) The things you need to clarify are:

1. What will the intervals in my GTF file represent?
2. What precise sequences will serve as the reference and provide interval coordinates for the GTF file?
3. How are those intervals to be determined?

It sounds like 3. is a big part of your question. If I'm reading things correctly, you want to annotate a particular class of genes on an influenza genome. Is that correct? Did you perform *de novo* sequencing of influenza viruses, and now hope to characterize variation in their HA genes? Stating what you hope to get out of your analysis would be helpful in finding a good way forward.

# Answer
> 1 votes
Since you know that your sequences in your fastas are your genes of interest, why not just dummy up a gtf where every sequence is a gene?

The point of a gtf is to tell software which genomic regions matter for a certain procedure, and which regions do not. But if you already know that your sequences are pared down to genes, then you know that the whole thing is what's relevant.

---
Tags: fasta, gtf, gff
---


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

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

---
Tags: database, bash
---


# Question
Title: Discordance in gene signature behavior between bulk and single-cell RNASeq
The objective of the following analysis is to identify an activation signature of a specific phenotype on bulk RNASeq and to apply it to single-cell RNA-Seq, in order to identify the population of activated cells in scRNA-Seq samples.

In the first step I used DESeq2 to identify genes that increase during different timepoints (I have 4 timepoints):

* Bulk samples are normalized using VST (Variance Stabilizing Transformation - DESeq2)
* For each timepoint, differentially expressed genes are computed with respect to the previous timepoint (i.e. 2 vs 1, 3 vs 2, 4 vs 3)
* Differentially expressed genes that increase with timepoint (log-FC \> 1, padj \< .05) are selected (n=24)

As a control, I plotted the activation score (computed with scanpy) based on the 24-gene signature applied to the bulk samples from which it has been extracted. As you can see, the score always increases from the first timepoint (blue) to the last (red).

However, when I use the same scoring procedure on a single-cell RNASeq sample containing the same populations, I think something may be wrong.

In the following the violin plot of the extracted activation score on an unsupervised clustering on single-cell data. Excluding cluster 4, which we identified as contamination, while cluster 1 and cluster 0 are supposed to have the correct degree of activation, clusters 2 and 3 have a lower degree of activation with respect to cluster 0 (which we identified as the first timepoint - lowest activation on bulk).

What may be the causes of this behavior? Right now the possible explanation I came out with is that:

Are there other possible explanations for this behavior, and how to check for them?

# Answer
> 3 votes
One explanation could be that your mapping of clusters to timepoints is not accurate. There are other methods you could look at for doing this, for example scMap, scPred, or Seurat v3 (disclosure: I am one of the Seurat developers).

Another possibility is that there is some problem in the calculation of the activation signature. You could instead look at the expression of the 24 genes themselves in each cluster (for example by drawing a heatmap).

The difference between bulk and single-cell here could also be biological. This could happen if the increase in expression of those genes in bulk was due to a change in cell type frequency over time rather than increased expression.

---
Tags: rna-seq, scrnaseq, normalization, tpm
---


# Question
Title: Mix globbing and wildcards when specifying rule input
Consider the following scenario.

```
├── sample-alice
│   ├── sequence_1.fastq
│   ├── sequence_2.fastq
│   ├── ...
│   └── sequence_n.fastq
└── sample-bob
    ├── sequence_1.fastq
    ├── sequence_2.fastq
    ├── ...
    └── sequence_m.fastq

```

I'm trying to write a single rule in my Snakefile that will preprocess each sample generically. The rule will execute a single command that takes all of Fastq files associated with the sample simultaneously. The samples do not have the same number of Fastq files.

If it was a single sample, I could do something like this.

```
rule preprocess:
    input: glob('sample-alice/*.fastq')
    output: 'alice-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

```

On the other hand, if there was a single Fastq file for each sample, I could use a wildcard to write a generic rule like this.

```
rule preprocess:
    input: 'sample-{samp}/sequence_1.fastq'
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

```

Is it possible to combine these two approaches to write a rule for handling multi-file samples generically?

# Answer
A nice solution comes from a colleague of mine, as seen in this Snakemake workflow. The trick is to access the wildcards programmatically using an anonymous (lambda) function. In my example above, it would be implemented as follows.

```
from glob import glob

rule preprocess:
    input: lambda wildcards: glob('sample-{samp}/*.fastq'.format(samp=wildcards.samp))
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

```
> 5 votes

---
Tags: snakemake
---


# Question
Title: bedtools: get name column of alignment file A
I'm using the `coverage` function of `bedtools` to check whether a certain set of regions (file B) has any overlap with known regions of epigenetic markers (file A). When that is the case, I want to access the name of that marker, which is given as the third column (name) in file A.

However, when I run `bedtools coverage -a fileA.bed -b fileB.bed`, I only get the columns of file B as an output, together with the overlap fraction. Is there a built-in way of accessing also the columns of the file A to which I aligned, or an easy way of cross-referencing to it?

# Answer
> 2 votes
I would recommend using bedtools intersect for that type of comparison. This will allow you to see the specific segments that overlap. bedtools coverage is used more for determining how many things cover a region but not which ones.

The documentation linked above will provide a much better desciption of intersect than I would be able to in this post but I would recommend using the `-loj` or `-u` depending on what you want.

```
-loj    Perform a “left outer join”. That is, for each feature in A report each overlap with B. If no overlaps are found, report a NULL feature for B.

-u  Write original A entry once if any overlaps found in B. In other words, just report the fact at least one overlap was found in B. Restricted by -f and -r.

```

# Answer
> 2 votes
If you're working with `fileB.bed` (a sorted BED file) and `fileA.txt` (a three-column text file, where the name is in the third column), then you could pipe single-base pseudo-BED elements to `bedmap`.

```
$ bedmap --echo --echo-map-id --delim '\t' fileB.bed <(awk -vFS="\t" -vOFS="\t" '{ print $1, $2, ($2+1), $3 }' fileA.txt) > answer.bed

```

The `--echo-map-id` option reports the list of marker IDs from `fileA.txt`, where the single-base marker interval from `fileA.txt` overlaps a region in `fileB.bed`.

If all the marker IDs are not unique (if there are duplicate marker IDs) it may make sense to instead use `--echo-map-id-uniq`.

This assumes that marker positions in the second column and that the intervals in `fileB.bed` are using a 0-based index scheme. If marker positions are 1-based and `fileB.bed` is 0-based, then a small adjustment is needed to make both inputs use the same indexing, to calculate correct overlaps:

```
$ bedmap --echo --echo-map-id --delim '\t' fileB.bed <(awk -vFS="\t" -vOFS="\t" '{ print $1, ($2-1), $2, $3 }' fileA.txt) > answer.bed

```

I mention this, as SNVs (such as those from VCF) can often be reported with a 1-based index scheme, which adds some nuance to doing set operations against BED files.

---
Tags: read-mapping, bedtools
---


# Question
Title: Alignment with arbitrary number of mismatches or gaps
I have 23bp long reads and want to find all possible alignments of them to the human genome (hg19, hg38) for an arbitrary number of mismatches (\<7), possibly also small indels. I've read in literature that people use bowtie2 for this, so I've tried the following:

```
echo "GGAGAACAAGTGGCAGATAGAGG" | bowtie2 --all -f -U- -x hg19 -r --no-hd --quiet --score-min C,-36 -p 3 -L 4

```

I'm not very experienced with bowtie2, and I hoped that this would give me all alignments of the piped-in sequence with up to 6 mismatches (6\*-6 = -36) to the hg19 genome which I previously indexed from fasta files. It however only gives me the correct full alignment:

```
0   0   chr2    173327389   255 23M *   0   0   GGAGAACAAGTGGCAGATAGAGG IIIIIIIIIIIIIIIIIIIIIII AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:23 YT:Z:UU

```

I've read that bowtie2 by design only allows a maximum of one mismatch per seed. Even if I implement the smallest possible seed length and interval (-L 4 -i C,1 -N 1), I only get this one result which I think is highly unlikely, given that I'm aligning to the whole genome.

The last three bp will have to match exactly - I was going to enforce this manually later on, maybe there's a better option for doing this. Taking the indels into account is also a secondary aim.

I'm fine with ~10k alignments per query for \<=6 mismatches plus the last 3bp matching exactly. Haven't done the combinatorics but that sounds realistic to me - question is how to make bowtie do it, or a different tool.

What would be the right way of achieving this using bowtie, or a different tool?

# Answer
This is not something mappers will be good at, because you're going to hit most of the genome. Bowtie won't do it, and I don't think LAST will either. The majority of fast mapping programs rely on seed-based mapping where the seed is essentially free of errors.

Doing only mismatches makes it much easier to determine, but it's still too relaxed. Consider that the example sequence you have provided is 10 mismatches away from `GGGGGGGGGGGGGGGGAGAGAGG`, and also 10 mismatches away from `AAAAAAAAAAAAAAAGAGAGAGG` (with the last 3 bp matching exactly)

Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:

```
for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors < 7)
            print gp

```
> 6 votes

# Answer
I've ended up using a package called BatMis for this, which can detect alignments with an arbitrary number of mismatches (not sure about indels). It is quite performant and yields about the order of magnitude of alignments that I was expecting (several thousands for alignment of a 20bp read to the whole human genome with up to 6 mismatches).
> 2 votes

# Answer
Did you see

> -k Report up to valid alignments per read or pair (default: 1).
> 0 votes

---
Tags: sequence-alignment, bowtie2
---


# Question
Title: I need some tips and suggestions for further analysis of NGS expression data (log2cpm)
I am a PhD student who inherited some log2cpm data of expression data from bulk kidney tissue from a UUO(unilateral urethral obstruction) experiment that tests a new drug. The sample material consists of:

* 6 x Ligated Kidney (Untreated
* 6 x Ligated Kideny (Treated)
* 3 x Unligated Kidney (Untreated)
* 3 x Unligated Kidney (Treated)

The previous study showed that this drug was effective against fribrosis, so the aim of my study is to investigate how this drug affects inflammation and mitochondrial function in the diseased kidney.

So I have filtered my original dataset based on MitoCarta v2 DB, so investigate how the drug affected mitochondrial function. I afterwards run a PCA analysis on log2cpm data which revealed the following:

The following figure shows PC1(84.51%) and PC2(4.69%). Red (Ligated, Untreated), purple (Ligated, Treated), green(healthy kidney, not treated), orange (Healthy kidney, treated).

I interpret this as there is a very small difference between treated and not treated ligated kidneys, however there is a small difference.

I have just started doing bioinformatics and I am very unsure which approach I should choose further, however I have though maybe do a highly variable genes analysis (HVG) see link, isolate the 500 most variable genes and run a heatmap with hieraichal clustering. I was also thinking of looking at the loadings in PC2, and take this further. The truth is that I am very unsure what is the right approach, and my research group does not have bioinformatician I can ask anymore.

I have the following questions:

* What kind of analysis you guys would suggest I use to investigate further (Hieraichal clustering, HVG etc?)
* Is there any alternative preprocessing I could use before I run the PCA, to make the analysis more accurate?
* Are there any packages that you can suggest?
* Any other tips?

Any other tips?

**I asked this same question on Biostars as well** https://www.biostars.org/p/367722/#367730

# Answer
Since you indicate already having performed differential expression and using gene set enrichment the next step is to look at what changes due to drug treatment in the pathways of interest while having a decent understanding of those pathways. Things like IPA or pathview (from Bioconductor) can aid in this, but at the end of the day you need to have a decent understanding of the pathways yourself to judge what seems to be like the best targets for future work.
> 1 votes

---
Tags: rna-seq, differential-expression, transcriptome, rna, gene-expression
---


# Question
Title: Error in formation of correlation matrix-(corrupt matrix -- dims not not match length)
I am working on a script to create a correlation matrix and mutual\_rank matrix from RNAseq gene counts for 57992 genes in 7027 samples. I have already tested this script for a smaller file (for 10,000 genes) and it is working, however, when I am running this script on my final input file, it is giving the following error:

```
 corrupt matrix -- dims not not match length (this error comes after gene to
 gene correlation step)

```

Script which i am running is:

```
E <- read.table("mdat_human_test.txt", header = TRUE, row.names=1)

# gene-to-gene correlation
C <- cor (t(E))
write.table(C, "human_correlation_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# correlation rank matrix
R <- t(apply(1-C, 1, rank))
write.table(R, "human_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# mutual rank matrix
M <- sqrt(R * t(R))
write.table(M, "human_mutual_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

q()

```

# Answer
> 3 votes
You can't compute a correlation matrix when you have more than ~46000 rows, since a standard R matrix can have a maximum of 2^31-1 values. Have a look at packages like bigcor. Alternatively, consider if having 3+ billion correlation values is really useful.

---
Tags: r, rna-seq, correlation
---


title: Getting data from fastq by generator

text:
# Question
Title: Getting data from fastq by generator
I have a task in a training that I have to read and filter the 'good' reads of big fastq files. I downsampled, got the code working, saving in a python dictionary. But turns out the original files are huge and I rewrite the code to give a generator. It did work for the down-sampled sample. But I was wondering if it's a good idea to get out all the data and filtering in a dictionary. Does anybody here a better idea? I am asking because I am learning python on my own and I am not sure if it makes sense.

I got some ideas from a code in Biostar:

```
import sys
import gzip

filename = sys.argv[1]

def parsing_fastq_files(filename):

    with gzip.open(filename, "rb") as infile:
        count_lines = 0
            for line in infile:
            line = line.decode()
            if count_lines % 4 == 0:
                ids = line[1:].strip()
                yield ids
            if count_lines == 1:
                reads = line.rstrip()
                yield reads
        count_lines += 1

total_reads = parsing_fastq_files(filename)
print(next(total_reads))
print(next(total_reads))

```

I now need to figure out to get the data filtered by using 'if value.endswith('expression'):' but if I use a dict for example, but that's my doubt because of the number of keys and values.

# Answer
> 2 votes
Have a function that returns a single-end read (or pair) and another that does any relevant filtering on that. Wrap both of them in another function that checks if you've hit the end of the file. That's both simple, will handle any size of file, and allows single and paired-end reads.

# Answer
> 2 votes
I do not have an example of fastqs with 'expression' as the last string in the header, so instead I will filter based off of a specific barcode used for demultiplexing. I have a gzipped fastq of 250 reads in size. Based on Devon's answer, I have first generated a fastq parser which stores only (header, seq) as a list of tuples, and then looks inside the list of tuples for a specific index. Some of the indexes vary by 1 nucleotide so those will be left out. I can prove this by printing the length of the 2 lists.

Here is the code:

```
import sys
import gzip

fastq = sys.argv[1]

def parser(filename):
    """
    Generate a list of tuples (header, read)
    """
    fastq_parsed = []
    try:

        with gzip.open(filename) as fq:
            header = next(fq)
            read = next(fq)
            fastq_parsed.append((header[1:-1], read[:-1])) # we don't want new-line chars or @ sign in reads
            while True:
                next(fq) # skip + line
                next(fq) # skip qscore
                header = next(fq) # store next read header
                read = next(fq) # store next read seq
                fastq_parsed.append((header[1:-1], read[:-1]))
    except:
        StopIteration # when we run out of reads, stop
    return fastq_parsed

def expression_finder(read_list):
    """
    Find headers that end with expression assuming the list format generated above.
    """
    filtered_reads = []
    for header, read in read_list:
        if (header.endswith(b'ATCACGAT')): #gzip files in bytes format, so we pass bytes object.
            filtered_reads.append((header.decode("utf-8"), read.decode("utf-8"))) # list items decoded to utf-8
    return filtered_reads

read_list = parser(fastq)
print(len(read_list))
subset_reads = expression_finder(read_list)
print(len(subset_reads))

```

Which returns:

```
(base) [dkennetz@nodecn203  test]$ python3.6 fastq_parser.py subset.fastq.gz
250
173

```

We can see that of the 250 reads, 173 have the correct index. I have hardcoded the filter in the second function, but you could also pass it as an argument to the function, like:

```
def expression_finder(read_list, filter):
    ...
        if header.endswith(filter):
            ...

subset_reads = expression_finder(read_list, b'expression')

```

which would also work. But include the `b''` in your string syntax because the first function still passes a bytes object.

---
Tags: python, fastq
---