# Question
Title: How can I improve the yield of MinION sequencing runs?
This is a frequently-asked question within the nanopore community. Oxford Nanopore currently claims that they are able to generate run yields of 10-15 gigabases (e.g. see here and here), and yet it's more common to see users only managing in the 1-5 gigabase range.

So why the big difference in yield?

# Answer
> 13 votes
I attended a talk by Josh Quick at PoreCampAU 2017, in which he discussed some common barriers to getting both good sequencing yield and long read length. It mostly boils down to being more careful with the sample preparation. Bear in mind that the MinION will still sequence a dirty sample, it will just be at a reduced yield. Here are my notes from that talk:

* The hardest thing about MinION sequencing is getting the sample in the first place
* There are lots of different sample types and extraction methods
* You can't get longer reads than what you put in; shit in leads to shit out
* DNA is very stable when not moving, but very sensitive to lateral damage
* The phenol chloroform method of DNA extraction is very good, and can be used with a phase-locked gel to make extraction easier
* A simple salt + alcohol extraction might be the best method for extraction (because it involves the least amount of work on the DNA)
* EDTA (e.g. as found in TE buffer, and many extraction kits) is not compatible with the rapid kit
* The most consistently good Nanopore runs produced by Josh's lab were 1D ligations runs on R9.4; the best overall run was a phenol-chloroform extraction + rapid kit
* John Tyson can tune himself out of a low-yield hole (via software)
* Getting small numbers of short reads is very important
* Suggested (and mostly untested) purification techniques: spin column (60-100kb); ethanol extraction (100-150kb), dialysis (150-250kb); low melting-point agarose plug (~1Mb, DNA extraction in situ)
* The nanopore protocol input is in nanograms, but should really be stated as molarity; the kit expects about 0.2 pmol input
* Picture molecules tethered to the surface of the membrane. You can then see that the flow cell density is independent of the sequence length
* Tapestation, Qubit and Nanodrop are all a good idea; a DNA sample that can pass all three tests will work well: no short-length shoulders by Tapestation, sufficient DNA by Qubit, high purity by Nanodrop
* RNA can interfere with sequencing; digesting RNA is recommended
* Freezing DNA is a really bad idea. The ice crystals are very good at chopping DNA up into small pieces
* DNA that is kept in the fridge is remarkably stable; can be kept for over two years (and probably indefinitely)

Update on yield:

We had a run in June 2018 which produced about 3 million reads from PCR-amplified mouse cDNA (read N50 of ~1kb), and that was in a situation where the it was the first sample prepared by someone, the input DNA amount was probably about 20% of what it should have been, and we suffered a dataloss incident. If that's what happened with a bad run, I have high hopes for our runs in the future.

Our next cDNA run in August 2018 produced about 7 million reads and a similar read N50 (i.e. total yield of about 7Gb). The run was interrupted by a Windows update during the first night (I had disabled automatic updates previously, but they had been re-enabled), and there have since been additional software and hardware updates since to improve sequencing yield. If we have a good flow cell and similar sample prep, I expect that the next cDNA run that we'll do should produce over 8 million reads, and possibly over 10 million.

# Answer
> 2 votes
My bet is that the nanopore folks have (of course!) done a ton of optimisation on two key things:

1. The DNA cleanup
2. The library preps

I think that for most DNA extractions cleanups can vary a lot, but something like a Blue Pippin will do an exceptional job if you have access to one and can use it (they won't work for SUPER long DNA fragments which can't move through a gel of course). The Blue Pippin also does your size selection, so you don't burn out pores sequencing lots of little stuff. This should help yield too.

Then it's all down to the library prep, and @gringer's answer has some good tips on that.

# Answer
> 1 votes
I am working with a fungus and am using CTAB method for DNA extraction and I used rapid kit for my initial library prep but it did not work and since than I have been using ligation kit and it has been giving me very good results. I managed to get data ranging from 4GB to 13GB.

---
Tags: fastq, minion
---


# Question
Title: How to align output of grep --color=always? (To QC fasta/fastq files)
Grepping out short sequences from a fasta or fastq file is a really useful way to look at sequencing data. Using the option `--color=always` makes this even more useful, as you can visualize where the sequences appear in sequencing reads. For example:

```
zgrep 'GATCGATC' --color=always My_Hiseq_Data.fq.gz

```

To make this more useful, how can I align the output of `grep --color=always` such that the beginning of the highlighted parts are aligned?

# Answer
Perhaps, `grep` is not the best tool to use in this case, but it should be in principle possible by using `grep` & `sed`. Here is an example showing three symbols around a match.

```
zcat My_Hiseq_Data.fq.gz | \
      grep -Eo '.{0,3}GATCGATC.*' | \
      sed -En 's/.*/   \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
      grep --color=always GATCGATC

```

Here are some explanations:

1. `zcat` will just decompress a file.
2. `grep` will output all the matches with 3-symbol context (or less). This is needed because a line can have multiple matches.
3. first `sed` command will add 3 symbols before all the matches.
4. second `sed` command will output the sequence with 3 surrounding characters on the left and max 3 characters on the right. This is an equivalent of `grep -Eo '.{3}GATCGATC.{0,3}'`.
5. Finally, now that all matches are aligned, `grep` will color your sequence.

This code is possible to adapt to include more symbols of context.

I would recommend doing this in another language, though. E.g., in Python:

```
#!/usr/bin/env python3
import sys
import re

HL = '\033[93m'
ENDC = '\033[0m'

filename = sys.argv[1]
seq = sys.argv[2]
context = 10

seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (context, sys.argv[2], context))

for line in open(filename):
    line = line.rstrip()
    for m in seqre.finditer(line):
        print(m.group(1).rjust(context),
              HL,
              m.group(2),
              ENDC,
              m.group(3),
              sep='')

```

# UPD

A more advanced Python version inspired by this answer with support for multiple files, gzip, etc.

```
#!/usr/bin/env python3
import sys
import re
import argparse
import gzip
from termcolor import  colored

def align_matches(f, seqre, context):
    # for every line
    for line in f:
        # get rid of newline character
        line = line.rstrip()
        for m in seqre.finditer(line):
            # for every match print right-aligned left-context,
            # match colored in red, and right-context
            print(m.group(1).rjust(context),
                  colored(m.group(2), 'red'),
                  m.group(3),
                  sep='')

# this function will check file extension to be able
# to open text files as well as gzip-compressed files
def myopen(fn):
    if fn.endswith('.gz'):
        return gzip.open(fn, 'rt')
    else:
        return open(fn, 'rt')

if __name__ == '__main__':
    # define script arguments
    parser = argparse.ArgumentParser(description='Aligned grep')
    parser.add_argument('pattern', type=str,
                        help='a valid regular expression')
    parser.add_argument('input', type=myopen, nargs='+',
                        help='files names, could be a gzip compressed file ending with .gz')
    parser.add_argument('--context', '-c', default=30, type=int,
                        help='number of context characters to show')
    args = parser.parse_args()

    # compile a regular expression with a match surrounded
    # by at most args.context number of characters
    seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (
        args.context, args.pattern, args.context))

    # for every input file run align_matches
    for f in args.input:
        align_matches(f, seqre, args.context)

```
> 6 votes

# Answer
I wrote a Perl script to do this. Just like `zgrep`, it can search through as many files as you give it. You can use the `-c` switch to tell it how much context (how many characters) you want around the match. So `-c10` will search for your target string and 10 characters on either side of it. If you don't use the `-c`, it will default to 155, which should include the entire sequence for most cases.

To use it, save as `alignString.pl` somewhere in your `$PATH` (`~/bin` should work if you're using Ubuntu, for example), make it executable with `chmod a+x ~/bin/alignString` and then run it like this:

```
alignString.pl GAAAATCAG file1.fastq.gz file2.fastq.gz file3.fastq.gz

```

For example:

The code:

```
#!/usr/bin/env perl

use warnings;
use strict;
use Term::ANSIColor;
use Getopt::Std;

my %opts;
getopts('c:',\%opts);
my $context = $opts{'c'} || 50;

## This subroutine does the alignment
sub align{
  ## Read the line and the target pattern
  my ($line, $pat) = @_;
  ## This will find 0 to $context characters on either side
  ## of the search pattern. Note that it will find the first
  ## occurrence of the pattern.
  $line =~ /(.{0,$context})$pat(.{0,$context})/;
  ## Save the strings to the left and right of the pattern
  my $left = $1;
  my $right = $2;
  ## If the left hand string is shorter than the value of $context,
  ## pad it with spaces so it can be aligned.
  while (length($left) < $context) {
    $left = " " . $left;
  }
  ## Return the padded left hand string, the pattenr colored red
  ## and the right hand string unchanged.
  return "$left" . colored($pat,"red") . "$right";
}

## The search pattern is the 1st argument
my $pattern = shift;
## Any remaining arguments are taken as target files to search
my @fastqFiles = @ARGV;

## Run the zgrep command
open(my $grepStream, '-|', "zgrep $pattern @fastqFiles");
## For each line returned by zgrep
while (my $line = <$grepStream>) {
  ## Align this line
  $line = align($line, $pattern);
  print "$line\n";
}

```
> 2 votes

---
Tags: sequence-alignment, fasta, fastq, illumina, bash
---


title: cDNA and alignment mapping

text:
# Question
Title: cDNA and alignment mapping
I am confused with RNA seq alignment. My understanding is that after Mature mRNA is isolated from the cell, it is then fragmented and using reverse transcriptase enzyme a cDNA copy is created which is further amplified. This cDNA is then mapped with the reference genome. If the cDNA does not have introns how does it map with the reference genome.

Can someone please help

# Answer
## tl;dr

Spliced aligners just make many small matches along the read pair or single read to determine splice junctions and output an alignment for the whole read, as long as it matches somewhere in the genome.

### mRNA and translation

In the cell, after transcription, introns are spliced out of RNAs to become mature RNA. The protein-coding sequence does have sequence that is chimeric from cis-located portions of the same chromosome that do not actually exist.

```
              |<<exon 1>>|  intron    |<<exon 2>|
Genome: ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

mRNA  : ----ABCDEFGHIJKLYZ123456789aaaaaaaaaaa

```

As you said, if we only even consider small stretches of DNA, the sequence `JKLYZ` does not exist in the genome, only in the mRNA. So, we need a special way to map mRNA-seq reads to the genome that correctly splits at splice junctions.

## How do we map mRNA-seq reads to a genome?

As a result, parts of the mRNA-seq reads will align well to one part of the genome, but only for part of the read. This requires a different alignment algorithm that considers multiple local high-score alignments for different parts of the read.

For example, see how a genomic aligner versus a splice aligner might treat the same mRNA-seq read on our example genome.

```
                         |<<exon 1>>|  intron    |<<exon 2>|
Genome:            ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

spliced aligner  : ------ABCDEFGHIJKL------------YZ123456789
genomic aligner  : ------ABCDEFGHIJKL

```

There is a possibility that the genomic aligner will not recognize that the second half of the read is an exon due to the large intron gap, and it may clip the alignment. Spliced aligners just take this into account and output several alignments for the same read regardless of gap size. The STAR algorithm, for example, uses many small matching pieces of a single read to determine where the splice junctions are. HiSat2 uses small approximate matches using a FM-index to find seeds for spliced alignment.
> 5 votes

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


# Question
Title: What i5 index should I use on the Illumina sample sheet for an unindexed p5 primer?
I have an upcoming run on a HiSeq X and most of the libraries in the pool have both i5 and i7 indices. However, some of the libraries were made with the IS4 p5 oligo and it is unindexed.

The IS4 indexing oligo is from the following paper. The sequence is below.

> Meyer, Matthias, and Martin Kircher. "Illumina sequencing library preparation for highly multiplexed target capture and sequencing." Cold Spring Harbor Protocols 2010.6 (2010): pdb-prot5448.

```
                                 miseq index
                              5' TCTTTCC 3'
                                 >>>>>>>
 5' AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTT 3'
    ^------------------^         ^--------------------^
    flowcell binding seq        adapter-binding seq

```

When we sequence these libraries dual-indexed on the MiSeq, the machine reads into the adapter and returns the index sequence `TCTTTCC`.

## Question

What will the index sequence be when the HiSeq reads it?

The more fundamental question is: How exactly do MiSeq and HiSeq machines read index sequences?

### side note

Here is what an index p5 and unindexed p5 molecule look like:

Update: Sept 2019 - After running into this problem again, an Illumina representative agreed with Devon Ryan and said that setting up a run for single indexing is the best option when there are unindexed p5 ends on one or more libraries.

# Answer
Illumina machines read the two indices in separate index sequencing reactions where they have sequencing barcodes specific for indices. If you have no index, then you should either not bother with creating either index read or mask that during demultiplexing. If you're mixing indexed samples with an unindexed sample, you can't really rely on any particular sequence being returned by the machine (i.e., just use the "undetermined indices" fastq file). The unindexed adapters tend to not have the full index sequencing primer binding site on them, so reads without index generally don't return much signal. In practice, they're likely to adopt barcode sequences of neighboring clusters. We see that a lot when we spike in PhiX. It has no barcode but you tend to get a small amount assigned to each sample either due to index hopping or (more likely) bleed through of signal from neighboring clusters (or small amounts of those clusters walking between nanowells and then their minority signal dominating during index sequencing).
> 1 votes

---
Tags: ngs, phylogenetics
---


title: Seurat for clustering bulk RNA-seq?

text:
# Question
Title: Seurat for clustering bulk RNA-seq?
Is it ever ok to use Seurat for clustering bulk samples?

I am looking at FPKM data from ~750 bulk RNA-seq samples generated using `Cufflinks`. As suggested for FPKM data, I manually input log transformed data to the `@data` slot \[`cd138_bm@data <- log(cd138_bm@raw.data + 1)`\] and skip the `NormalizeData()` function. I then use functions `FindVariableGenes`, `ScaleData`, `RunPCA`, `FindClusters`, `RunTSNE`, `FindAllMarkers` in their usual ways to find clusters & cluster markers. My clustering results are quite reasonable, and reflect published work clustering similar samples.

What are the potential pitfalls of using these Seurat functions on bulk data? In `FindAllMarkers`, would you recommend I use the "negbinom" test? (currently using wilcox) Any other arguments you would recommend changing from the default?

# Answer
> 2 votes
I'm not sure Seurat is the best tool for this as it was developed for single cell RNA seq data and there are a few intricacies of that type of data that are very different from bulk RNA seq. For bulk there are really good packages available and corresponding workflows, e.g. limma, edgeR and DESeq2.

The main problems with Seurat for bulk RNA-seq:

1. Seurat expects counts as input - FPKM are not counts nor are they log counts or log norm counts. It's a specific normalisation method that takes into account gene length and library size and breaks the link between gene counts and variance.
2. Because of the distribution of reads from a single cell that you get from the sequencer Seurat scales all unique reads (i.e. UMIs) to 10,000 per cell and usually regresses out the nUMI count. Again, becuase you've got FPKMs this makes no sense.
3. The other functions you can use on bulk data, i.e. tSNE, finding clusters and markers (which is essentially DGE analysis). The one thing to note for the latter is that again because you've got FPKMs they don't follow the negative binomial distribution. Your best bet is limma trend but I'm not sure it's included in Seurat.

# Answer
> 1 votes
FYI - heard from the author of Seurat, see this github issue

His suggestion to use DESeq won't work in my case, because DESeq requires counts data which I don't have. But it is the green light to generally use the package on bulk data. I am using test.use = "wilcox" in FindMarkers() because it's non-parametric...

---
Tags: rna-seq, seurat, clustering
---


title: 3D PCA group labelling

text:
# Question
Title: 3D PCA group labelling
I would like to make a 3D PCA but not sure how to label group wise which i can do for 2D PCA

```
library(ggplot2)
test<- read.csv("NEW_RBP/TF_DISEASE.txt",header = T,row.names = 1,sep = '\t')

Sample = c(rep("HSC",4),rep("Blast",11),rep("LSC",8))
Stage=c(rep("HSC",4),rep("Blast",11),rep("LSC",8))

df_final=data.frame(Sample,Stage)
df_final

pca_data=prcomp(t(test), center=TRUE, scale=TRUE)
pca_data_perc=round(100*pca_data$sdev^2/sum(pca_data$sdev^2),1)

df_pca_data=data.frame(PC1 = pca_data$x[,1],
                       PC2 = pca_data$x[,2], sample = colnames(test))
myColors <- c("red", "black", "blue2")
sample <- as.vector(df_pca_data$sample)
ggplot(df_pca_data, aes(PC1,PC2, colour = factor(Stage)))+#,label=rownames(t(test)))+
  geom_point(size=30)+
  geom_text(aes(label=sample),color="white",size=6,angle =0,parse = TRUE,  family="Bookman", fontface="bold")+
  labs(x=paste0("PC1 (",pca_data_perc[1],")"), y=paste0("PC2 (",pca_data_perc[2],")"))+
  theme_minimal(base_size=20) +
  theme(axis.text.x=element_text(size=rel(3), angle=90))+
  theme(axis.text.y=element_text(size=rel(3), angle=90))+
  theme(axis.title.x = element_text(colour="grey20",size=55,angle=0,hjust=.5,vjust=0,face="plain"))+
  theme(axis.title.y = element_text(colour="grey20",size=55,angle=90,hjust=.5,vjust=.5,face="plain"))+
  scale_color_manual(values=myColors)

```

The above code give me the figure as i have attached

For 3D PCA I tried this

```
scores = as.data.frame(pca_data$x)
head(scores)

plot3d(scores[,1:3], col=c(1:4), size=20, type='p',
       xlim = c(-50,50), ylim=c(-50,50), zlim=c(-50,50))
text3d(scores[,1]+2, scores[,2]+1, scores[,3]+1,
       texts=c(rownames(scores)), cex= 0.7, pos=3)

```

I get something like this as in case of 2D PCA i can label lets say all of the HSC sample as black and so on..Im not sure how to do that for 3D as I tried using factor i couldn;t do

Any help or suggestion would be highly appreciated

# Answer
> 5 votes
You can easily color 3D pca plots in R based on the code given below:

```
library("scatterplot3d")
colors <- c("#999999", "#E69F00", "#56B4E9") # Number of color according to the number of groups
colors <- colors[as.numeric(iris$Species)] # you can put here the column containing the name of population or sample etc.
pca1 <- prcomp(iris[, -5])
s3d <-scatterplot3d(pca1$x[, 1], pca1$x[, 2],pca1$x[, 3],xlab="Sepal.length",ylab="Sepal.width", zlab="Petal.length", pch = 16,color=colors)
legend("right", legend = levels(iris$Species),
  col =  c("#999999", "#E69F00", "#56B4E9"), pch = 16)

```

The generated graph is given below:

For further details about how to manipulate plots in 3D have a look at this site.

**Edited Response** Using the data given below I plotted a 3D pca that may help solve your problem,

```
"Col1" "Col2" "Col3" "Col4" "colend"
"H1" 5.1 3.5 1.4 0.2 "HSC"
"H2" 4.9 3 1.4 0.2 "HSC"
"H3" 4.7 3.2 1.3 0.2 "HSC"
"H4" 4.6 3.1 1.5 0.2 "HSC"
"H5" 5 3.6 1.4 0.2 "HSC"
"H6" 5.4 3.9 1.7 0.4 "HSC"
"H7" 4.6 3.4 1.4 0.3 "HSC"
"H8" 5 3.4 1.5 0.2 "HSC"
"H9" 4.4 2.9 1.4 0.2 "HSC"
"H10" 4.9 3.1 1.5 0.1 "HSC"
"H11" 5.4 3.7 1.5 0.2 "HSC"
"H12" 4.8 3.4 1.6 0.2 "HSC"
"H13" 4.8 3 1.4 0.1 "HSC"
"H14" 4.3 3 1.1 0.1 "HSC"
"H15" 5.8 4 1.2 0.2 "HSC"
"B1" 5.7 4.4 1.5 0.4 "blast"
"B2" 5.4 3.9 1.3 0.4 "blast"
"B3" 5.1 3.5 1.4 0.3 "blast"
"B4" 5.7 3.8 1.7 0.3 "blast"
"B5" 5.1 3.8 1.5 0.3 "blast"
"B6" 5.4 3.4 1.7 0.2 "blast"
"B7" 5.1 3.7 1.5 0.4 "blast"
"B8" 4.6 3.6 1 0.2 "blast"
"B9" 5.1 3.3 1.7 0.5 "blast"
"B10" 4.8 3.4 1.9 0.2 "blast"
"B11" 5 3 1.6 0.2 "blast"
"B12" 5 3.4 1.6 0.4 "blast"
"B13" 5.2 3.5 1.5 0.2 "blast"
"B14" 5.2 3.4 1.4 0.2 "blast"
"B15" 4.7 3.2 1.6 0.2 "blast"
"B16" 4.8 3.1 1.6 0.2 "blast"
"B17" 5.4 3.4 1.5 0.4 "blast"
"B18" 5.2 4.1 1.5 0.1 "blast"
"B19" 5.5 4.2 1.4 0.2 "blast"
"B20" 4.9 3.1 1.5 0.2 "blast"
"L1" 5 3.2 1.2 0.2 "LSC"
"L2" 5.5 3.5 1.3 0.2 "LSC"
"L3" 4.9 3.6 1.4 0.1 "LSC"
"L4" 4.4 3 1.3 0.2 "LSC"
"L5" 5.1 3.4 1.5 0.2 "LSC"
"L6" 5 3.5 1.3 0.3 "LSC"
"L7" 4.5 2.3 1.3 0.3 "LSC"
"L8" 4.4 3.2 1.3 0.2 "LSC"
"L9" 5 3.5 1.6 0.6 "LSC"
"L10" 5.1 3.8 1.9 0.4 "LSC"
"L11" 4.8 3 1.4 0.3 "LSC"
"L12" 5.1 3.8 1.6 0.2 "LSC"
"L13" 4.6 3.2 1.4 0.2 "LSC"
"L14" 5.3 3.7 1.5 0.2 "LSC"
"L15" 5 3.3 1.4 0.2 "LSC"

```

The code is as follows:

```
library("scatterplot3d")
colors <- c("#999999", "#E69F00", "#56B4E9") # Number of color according to the number of groups
colors <- colors[as.numeric(data.c1$colend)] # you can put here the column containing the name of population or sample etc.
pca1 <- prcomp(data.c1[, -5]) # PCA on columns except the last column
s3d<-scatterplot3d(pca1$x[, 1], pca1$x[, 2],pca1$x[, 3],grid=TRUE,xlab="PC1",ylab="PC2", zlab="PC3", pch = 16,color=colors)
legend("right",legend = levels(data.c1$colend),col =  c("#999999", "#E69F00", "#56B4E9"), pch = 16,inset=-0.04,bty="n")
text(s3d$xyz.convert(pca1$x[, 1:3]), labels = rownames(data.c1),cex= 0.7, col = "black",pos=2.5)

```

In this way we get following plot:

Hope this helps!!!

---
Tags: r, statistics, visualization, pca
---


# Question
Title: How to get the corrected matrix after SVA batch effect correction
I ran `SVA` to remove batch effects for my `bulk RNAseq` experiments, but now I need to somehow correct my data matrix in order to run `pca`, `mds`. I am using `DESeq2` for the analysis. Here is the code that I got now:

```
dds <- estimateSizeFactors(dds)
dat = counts(dds, normalized = TRUE)
idx = rowMeans(dat) > 1
dat = dat[idx,]
mod = model.matrix(~Group, colData(dds))
mod0 = model.matrix(~1, colData(dds))

nsv = num.sv(dat, mod)
svseq = svaseq(dat, mod, mod0, n.sv = nsv)

dds_sva = dds
dds_sva$SV1 = svseq$sv[,1]
design(dds_sva) = ~ SV1 + Group
dds_sva = DESeq(dds_sva)

```

So, I am wondering now, how to use `svseq` object to plot, say, the corrected version of `MDS plot`.

# Answer
> 2 votes
If you want to plot the "corrected" expression, you will need to remove the variation introduced by these surrogate variables. Removing the expression affected can introduce some bias too and it is usually not recommended (despite comBat doing so). You should apply linear algebra, you can look at here is an example how to do it:

```
library("corpcor")
s <- corpcor::fast.svd(t(scale(t(dat), center = TRUE, scale = TRUE)))
pcSds <- s$d
pcSds[1] <- 0
svdexp <- s$u %*% diag(pcSds) %*% t(s$v)
colnames(svdexp) <- colnames(dat)

```

But if your surrogate variables do have a big enough influence on your MDS plot to be observable (IMHO) there would be too much batch effect to work with the data.

---
Tags: r, rna-seq, deseq2, batch-effects, sva
---


title: RSVSim insertions from chr1 to chr2

text:
# Question
Title: RSVSim insertions from chr1 to chr2
# edit explaining python tag

I would still rather have a solution based on the RSVSim package in R, but while waiting for someone who might have an answer I wanted to look at other solutions as well.

If a result like what metadata(sim) produces can be achieved in python as well, I can do it without RSVSim. I need especially the numbers given there, so I do NOT need "Copied BpSeqA BpSeqB\_5prime BpSeqB\_3prime". The numbers need to reflect the original position of insert into chr2, because I need to reflect this back to the original (unedited) chr2.

# end edit explaining pyton tag

I tried using the BioConductor forums for asking this question Question on BioConductor. However I do not seem to get an answer there and after some help with other tools on this website, i thought I could also ask the question here, in the hope that someone knows something about it.

I have a reference genome, and I need to put in some structural variations. It is really easy to do in RSVSim, and used it quite a lot already.

However, when processing INS, I found out that it just copies a DNA seq from one part, and puts it somewhere else in the same DNA seq (Bascially a DUP then), or if you specify another an extra DNAseq, to cut from, It randomly cuts from either chr1 or chr2 and puts it in the other. This creates deletions in the DNA seq I want to investigate, which is not what I want.

Does anyone know how to get this kind of result?

```
width seq                                                                                names
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT    chr1
[2]    40 GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC    chr2

```

to this:

```
width seq                                                                                names
[1]    32 AAAAAAAAAAAAAAAATTTTTTTTTTTTTTTT                    chr1
[2]    48 GGGGGGGGGGGGGGAAAATTTTGGGGGGCCCCCCCCCCCCCCCCCCCC    chr2

```

And do it 30 times, randomly in the genome? (so take a seq from chr1 and put an insertion in chr2.

I tried this and more already for days, and I just cannot figure it out.

```
    library(RSVSim)

    seq_random <- readDNAStringSet("/path/random.fasta", "fasta")
    seq_ref <- readDNAStringSet("/path/reference1MB.fasta", "fasta")
    genome2 = DNAStringSet(c(seq_random, seq_ref))
    names(genome2) = c("chr1","chr2")

    y <- c(30, 30, 30, 30, 30, 50, 50, 50, 50, 50, 100, 100, 100, 100, 100, 500, 500, 500, 500, 500, 1000, 1000, 1000, 1000, 1000, 5000, 5000, 5000, 5000, 5000)

    length_seq = width(genome2[2])
    knownInsertion = GRanges(IRanges(0,length_seq), seqnames="chr2")
    knownInsertion = GRanges(IRanges(0,length_seq), seqnames="chr1", chrB="chr2", chrA="chr1")
    knownInsertion = GRanges(seqnames="chr1", chrB="chr2")
    sim = simulateSV(output='/folder/of/output/', genome=genome2,
                     ins = 30, sizeIns=y, regionsIns=knownInsertion)

```

## EDIT for comment section

The simple genome example:

```
genome = DNAStringSet(c("AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT","GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))
names(genome) = c("chr1","chr2")

```

I believe I already tried those parameters. For example:

```
knownInsertion = GRanges(IRanges(16,25), seqnames="chr1", chrB="chr2")
sim = simulateSV(output=NA, genome=genome, ins = 1, regionsIns=knownInsertion, bpSeqSize=6, random=TRUE)

```

results in:

```
> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT                                                                                chr1
[2]    40 GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC                                                                                chr2
> metadata(sim)
$insertions
         Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1      3   12 chr1      6   15   10  FALSE  AAAAA        AAAAAA        AAAAAA

```

so nothing happened at all. If I keep it simple, you get the result I mentioned before, where it is sometimes chr1-\>chr2 and sometimes chr2\<-chr1

```
sim = simulateSV(output=NA, genome=genome, ins = 3, sizeIns = 5, bpSeqSize=6, random=TRUE)

> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    40 AAAAAACCCCCAAAAATTTTTTTTTTTTTTAAAATTTTTT                                                                                chr1
[2]    40 GGGGGGGGAAAAAGGGGGGGGGGGGCCCCCCCCCCCCCCC                                                                                chr2
> metadata(sim)
$insertions
         Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1     17   21 chr1     36   40    5  FALSE AAATTT        TTTAAA        AATTTT
2 insertion_2 chr1      1    5 chr2      9   13    5  FALSE               GGGAAA        AAAGGG
3 insertion_3 chr2     33   37 chr1     12   16    5  FALSE CCCCCC        AAACCC        CCCAAA

```

In this case I dont even know what happened exactly, because an uneven number of Ins, should make one of the chromosomes shorter or longer?

Finally, when setting copy on (so no cut/paste)

```
sim = simulateSV(output=NA, genome=genome, ins = 3, sizeIns = 5, bpSeqSize=6, random=TRUE, percCopiedIns=1.00)
sim
metadata(sim)

> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    50 AAAATTTTTAAAAAAAAAAAAAAAATTTTTTATTTTTTTTTTTTTTTTTT                                                                            chr1
[2]    45 GGGGGGGGGGGGGGGGGAAAAAGGGCCCCCCCCCCCCCCCCCCCC                                                                                chr2
> metadata(sim)
$insertions
         Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1 insertion_1 chr1     36   40 chr1      5    9    5   TRUE               AAATTT        TTTAAA
2 insertion_2 chr1     20   24 chr1     27   31    5   TRUE               TTTATT        TTTTTT
3 insertion_3 chr1     13   17 chr2     18   22    5   TRUE               GGGAAA        AAAGGG

```

Again. all the examples above have random chrA and chrB. Setting these specifically, only seems possible with the GRanges and IRanges, for which I tried tons of versions already.

It should be a very simple thing according to the manual, but actually doing it doesn't seem to get the desired results, while my version of Insertions is just the normal one in my opinion? (why would you want to create extra DELS or DUPS while making insertions?).

On their own support page, it is only possible to ask a question on the forums, but they never answer it (only look at it)

# Answer
Answer eventually found on the BioConductor post: https://support.bioconductor.org/p/114354/#114585

It basically comes down to not using RSVSim, but other R packages.

Final code snippit:

```
genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

sampleFromGenome <- function(genome, wd) {
  chrProb <- width(genome) / sum(width(genome))
  chr <- sample(names(genome), n, TRUE, prob = chrProb)
  start <- runif(n, 1, nchar(genome[chr]) - wd)
  GRanges(IRanges(start, width = wd), seqnames = chr)
}

n <- 2
wd <- sample(5:10, n, TRUE)
insertFrom <- sampleFromGenome(genome[1], wd)
insertAt <- sampleFromGenome(genome[2], 0)

#writeTable(insertFrom, stdout()) . # replace stdout with a file name

sequence <- genome[insertFrom]
at <- split(insertAt, seqnames(insertAt))[names(genome[2])]
atSequence <- split(sequence, seqnames(insertAt))[names(genome[2])]

replaceAt(genome[2], unname(ranges(at)), atSequence)

  A DNAStringSet instance of length 1
    width seq                                                                                names
[1]    54 GGGGGGGGGGGGGGGGGGTTTTTTTTTGGCCCCCCCCCCCCTTTTTCCCCCCCC                                                                        chr2

```

# EDIT with RSVSim solution!

Finally found a solution based on RSVSim

```
genome = DNAStringSet(
  c(chr1 = "AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT",
    chr2 = "GGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCC"))

length_seq_chr1 = width(genome[1]) - 10
length_seq_chr2 = width(genome[2]) - 10
sizes <- c(3,3,5,5)
start_chr1 <- c("list", length(sizes))
start_chr2 <- c("list", length(sizes))
end_chr1 <- c("list", length(sizes))
index = 1
for(i in sizes){
  newstartchr1 <- sample(1:length_seq_chr1,1)
  newstartchr2 <- sample(1:length_seq_chr2,1)
  start_chr1[[index]] <- newstartchr1
  start_chr2[[index]] <- newstartchr2
  end_chr1[[index]] <- newstartchr1 + (i-1)
  index <- index + 1}
start_chr1 <- as.integer(start_chr1)
end_chr1 <- as.integer(end_chr1)
start_chr2 <- as.integer(start_chr2)

knownInsertion2 = GRanges(IRanges(c(start_chr1),c(end_chr1)),seqnames="chr1", chrB="chr2", startB=c(start_chr2), copied=TRUE)
knownInsertion2

sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
sim
metadata(sim)

```

Gives as an result:

```
> knownInsertion2
GRanges object with 4 ranges and 3 metadata columns:
      seqnames    ranges strand |        chrB    startB    copied
         <Rle> <IRanges>  <Rle> | <character> <integer> <logical>
  [1]     chr1  [21, 24]      * |        chr2        30      TRUE
  [2]     chr1  [16, 19]      * |        chr2        21      TRUE
  [3]     chr1  [ 6, 11]      * |        chr2        10      TRUE
  [4]     chr1  [ 9, 14]      * |        chr2        27      TRUE
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
>
> sim = simulateSV(output=NA, genome=genome, regionsIns=knownInsertion2, bpSeqSize=5, random=FALSE, verbose=FALSE)
> sim
  A DNAStringSet instance of length 2
    width seq                                                                                names
[1]    40 AAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTTTT                                                                                chr1
[2]    60 GGGGGGGGGAAAAAAGGGGGGGGGGGAAAACCCCCCAAAAAACCCTTTTCCCCCCCCCCC                                                                  chr2
> metadata(sim)
$insertions
  Name ChrA StartA EndA ChrB StartB EndB Size Copied BpSeqA BpSeqB_5prime BpSeqB_3prime
1    1 chr1     21   24 chr2     30   33    4   TRUE                 CCTT          TTCC
2    2 chr1     16   19 chr2     21   24    4   TRUE                 GGAA          AACC
3    3 chr1      6   11 chr2     10   15    6   TRUE                 GGAA          AAGG
4    4 chr1      9   14 chr2     27   32    6   TRUE                 CCAA          AACC

```

To make this customizable change:

sizes \<- c(whatever, numbers, you, like)

output=NA (to a file path prefix to get a .csv file with the metadata)

genome = DNAStringSet() To replace the sample sequence

bpSeqSize = number (to change how many of the flanking region is shown in the metadata)
> 2 votes

---
Tags: r, python, bioconductor, structural-variation, simulated-data
---


# Question
Title: API/website for blasting peptide sequence against database of all plant proteomes
I would like to blast a peptide sequence against all available plant proteome databases, i.e. a `blastp` metaproteome analysis.

What tools are available? Which database can I use?

# Answer
> 1 votes
A simple BLAST on the UniProt website does the trick (target database *plants*).

Best matches are subsequences `ETVEREKE` and `EKKRRETV` in *Green alga* and *Eelgrass* respectively.

---
Tags: proteins, blast, metagenome
---


# Question
Title: Total reads aligning to each reference within a bam file
I have two PCR amplicons that have been multiplexed and sequenced using the nanopore minion.

I have aligned the fastq reads using minimap2 with a reference file containing both amplicon sequences and generated a bam file that I have viewed using IGV.

I am looking for a way to generate some simple summary statistics.

In particular, is there a way to extract the total number of fastq reads aligning to each amplicon reference from the bam file?

# Answer
> 8 votes
The quick way to get the number of alignments on each reference is

```
samtools idxstats my_bam.bam

```

Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of *alignments* per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us this as generally I'm after a rough approximation, rather than an accurate number.

In theory, only one alignment for each read should be marked as primary, so the following should give you what you need quickly and at low memory usage:

```
samtools view -bF 2304 my_bam.bam > primary_only.bam
samtools index primary_only.bam
samtools idxstats primary_only.bam

```

# Answer
> 2 votes
This one-liner below will work better for long reads than `samtools flagstat` in that it only counts the primary alignment for each read and `samtools flagstat` doesn't seem to calculate some stats for long reads. I have never seen `samtools flagstat` output stats on a per-reference basis, but am curious if so!

This answer filters out secondary and supplementary alignments for your reads (`-F 2304`) that have some alignment to both amplicon reference and just keep the best one. This might give a more accurate idea of how many reads of each amplicon are in the library.

```
samtools view -F 2304 myfile.bam | awk -F $'\t' '{a[$1, $3]++} END{for (i in a) {split (i, sep, SUBSEP); print sep[1], sep[2], a[i]}}' | uniq | awk '{print($2)}' | uniq -c | sort -k1 -nr

```

adapted from this.

# Answer
> 1 votes
This can be done using

```
   samtools flagstat your_bam_file

```

---
Tags: sequence-alignment, sam, nanopore, minion
---


# Question
Title: Changing a wide range of colours to a limited gradient
I returned a FeaturePlot from Seurat to ggplot. My plot has a weird range of colours as below

I produced this plot by this code

```
> head(mat[1:4,1:4])
             s1.1        s1.2 s1.3       s1.4
DDB_G0267178    0 0.009263254    0 0.01286397
DDB_G0267180    0 0.000000000    0 0.00000000
DDB_G0267182    0 0.000000000    0 0.03810585
DDB_G0267184    0 0.000000000    0 0.00000000
>

I have converted expression matrix to a binary matrix by 2 as a threshold

mat[mat < 2] <- 0
mat[mat > 2] <- 1

> head(exp[1:4,1:4])
             s1.1 s1.2 s1.3 s1.4
DDB_G0267382    0    0    0    1
DDB_G0267438    0    0    0    1
DDB_G0267466    0    0    0    0
DDB_G0267476    0    0    1    0
>
> exp=colSums(exp)
> exp=as.matrix(exp)
> colnames(exp)="value"
> exp=as.data.frame(exp)
> cc <- AddMetaData(object = seurat, metadata = exp)
> cc <- SetAllIdent(object = cc, id = "value")
> TSNEPlot(object = cc, do.return= TRUE)

```

How I can convert this range to a gradient of colours for example in 8-18, 18-28, 28-38, 38-48 range in blue to yellow please? Something like below

Thank you for any help

Then by ggplot now I scaled my colours but I don't like my clusters as so and I don't know how to retain my clusters as a featureplot by this new color gradient

```
> head(cc@meta.data)
     nGene    nUMI    orig.ident res.0.7 CELL STAGE GENO dataset stage.nice celltype value
s1.1  4331  373762 SeuratProject       0 s1.1   H16   WT       1        H16        0    34
s1.2  5603 1074639 SeuratProject       0 s1.2   H16   WT       1        H16        0    26
s1.3  2064   49544 SeuratProject       0 s1.3   H16   WT       1        H16        0    27
s1.4  4680  772399 SeuratProject       1 s1.4   H16   WT       1        H16        1    29
s1.5  3876  272356 SeuratProject       1 s1.5   H16   WT       1        H16        1    21
s1.6  2557  122314 SeuratProject       0 s1.6   H16   WT       1        H16        0    31

> ggplot(as.data.frame(cc@meta.data), aes(x = cc@meta.data$CELL, y = cc@meta.data$res.0.7, colour =cc@meta.data$value)) +
+     geom_point(size = 5) +
+     scale_colour_gradient(low = "yellow", high = "blue")

```

By below code I obtained a tsne in link

```
> cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))
>
> TSNEPlot(cc, colors.use=cols)

```

enter link description here

Now I want to know how I could convert this range to a 8-18, 18-28, 28-38, 38-48 colour range as a gradient of blue to yellow?

# Answer
> 3 votes
If you would like to color discrete intervals on a gradient as opposed to having a continuous gradient (like your second plot), use this approach.

It is similar to the approach in the answer I posted with the continuous scale, but we simply break up the continuous scale in to intervals and color them by these intervals.

```
#generate values for testing purposes, one value for each cell, add to object
value <- sample(seq(from=8, to=46, by=1), size = length(rownames(exp@meta.data)), replace=TRUE)
exp <- AddMetaData(object=exp, metadata=value, col.name="value")

#encode the continuous values as factors, determined by the interval they fall in to
value_breaks <- cut(exp@meta.data$value, breaks = c(8,18,28,38,46), include.lowest=TRUE, right=FALSE)

#name the breaks by cell so they can be added by AddMetaData(), add them
names(value_breaks) <- rownames(exp@meta.data)
exp <- AddMetaData(object=exp, metadata=value_breaks, col.name="value_breaks")

exp <-  SetIdent(exp, ident.use=exp@meta.data$value_breaks)

#create a gradient between specified intervals
#length.out is the number of different colors, number of factor levels
cols <- scales::seq_gradient_pal(low="blue", high="yellow", space="Lab")(seq(from=0, to=1,length.out=4))

TSNEPlot(exp, colors.use=cols)

```

The label in the legend are classic mathematical notation for intervals. You can add your own labels if you wish in the `cut()` function.

# Answer
> 3 votes
**TSNEPlot()**

`TSNEPlot()` will always treat your variables as discrete. My approach is to manually generate a gradient with unique colors for each factor level and pass it to the `cols.use` argument in `TSNEPlot()`.

```
#generate values for testing purposes, one value for each cell
value <- sample(seq(from=8, to=48, by=1), size = length(rownames(unfiltered_cca@meta.data)), replace=TRUE)
names(value) <- rownames(exp@meta.data)
exp <- AddMetaData(object=exp, metadata=value, col.name="value")
exp <-  SetIdent(exp, ident.use=exp@meta.data$value)
#create a gradient between specified colors, multiply by sequence to get appropriate length
#length.out is the number of different colors, number of factor levels
cols <-  scales::seq_gradient_pal(low="beige", high="red", space="Lab")(seq(from=0, to=1,length.out=48))

TSNEPlot(exp, colors.use=cols)

```

I would present this with no legend:

```
TSNEPlot(exp, colors.use=cols, no.legend = TRUE)

```

**FeaturePlot()**

You can also simply use `FeaturePlot()` instead of `TSNEPlot()` to visualize the gradient. Using the same data as above:

```
FeaturePlot(object = exp, features.plot = "value", reduction.use = "tsne", no.legend = FALSE, cols.use = c("beige", "red"))

```

You ask for a continuous scale, but this is not what is shown in your second plot. You have a different color for each discrete range of values. I assume here that you want a continuous gradient. To achieve something like the second plot you provided requires a different approach.

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


# Question
Title: Fast filtering of intervals not falling within a certain distance from known genes
I would like to filter a bed file with intervals, ie. in the format of:

```
chr1    13800   14301
chr1    15500   16001
chr1    19400   19901
chr1    22800   23301

```

In particular, I want to filter out intervals that fall far from known genes, let's say \> 1KB.

I was thinking of querying each interval (in parallel) and check through a known set of coordinates from a gtf, e.g. from `Homo_sapiens.GRCh38.93.gtf` by filtering for gene coordinates:

```
1       havana  gene    168317497       168376876       .       +       .       gene_id "ENSG00000227722"; gene_version "2"; gene_name "AL022100.1"; gene_source "havana"; gene_biotype "transcribed_processed_pseudogene";
1       mirbase gene    168375524       168375621       .       +       .       gene_id "ENSG00000207974"; gene_version "1"; gene_name "MIR557"; gene_source "mirbase"; gene_biotype "miRNA";
1       havana  gene    168400829       168495685       .       -       .       gene_id "ENSG00000228697"; gene_version "2"; gene_name "AL023755.1"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168401483       168407274       .       -       .       gene_id "ENSG00000283255"; gene_version "1"; gene_name "AL022100.2"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168449672       168451211       .       -       .       gene_id "ENSG00000215833"; gene_version "3"; gene_name "QRSL1P1"; gene_source "havana"; gene_biotype "processed_pseudogene";

```

However, since the bedfile is quite big I am looking for the fastest solution since checking each interval may be quite expensive.

# Answer
> 5 votes
Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.

1. Convert GTF annotations to a sorted BED file of genes:

```
$ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed

```
2. Sort your intervals, if unsorted:

```
$ sort-bed intervals.unsorted.bed > intervals.bed

```
3. Filter intervals against 1kb-padded genes:

```
$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - > intervals.filtered.bed

```

We first `--merge` genes to simplify dealing with overlapping genes and further speed downstream filtering.

The file `intervals.filtered.bed` will contain intervals from `intervals.bed` which lie outside of a 1kb-window around merged gene space.

---

To run all of this as a one-liner (assuming use of the `bash` shell):

```
$ bedops --merge <(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 <(sort-bed intervals.unsorted.bed) - > intervals.filtered.bed

```

This creates the same output result, but avoids creation of intermediate files and so should run even faster.

---

**Note**: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using `--range` takes care of ranges that would have a negative start position, setting their start coordinates to `0`. Bounds aren't really relevant here, in any case, because the `bedops --not-element-of` step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.

# Answer
> 5 votes
In one line, using `bedtools`

```
zcat  Homo_sapiens.GRCh38.93.gtf.gz \
 | awk '$3=="gene"' \
 | bedtools slop -b 10000 -g contigs.tsv -i - \
 | bedtools intersect -u -a intervals.bed -b -

```

This first takes the genes and filters them for the third column being gene. This is important because all `GTF` lines contain the word `gene` as `gene_id` is a manditory attribute for `GTF`.

Next it passes it to slop, which adds 1 kb either side (ensuring not to go below 0 or off the end of the chr: this isn't that important because we are never outputting these regions anyway, but `bedtools` requires it). We tell `bedtools` to work on the standard in by using `-i -`.

Finally we intersect these gene regions with the `intervals.bed`. `-u` says output the line from `-a` if it overlaps with at least one entry in `-b` and passing `-b -` means read file b from the standard in.

This will only work if your `GTF` contains `gene` lines, which the GENCODE/Ensembl ones do. However, it is interesting to note that these GTFs are not actaully valid GTF files, as to be valid a GTF file must have both a `gene_id` and `transcript_id` attribute, and the gene lines in the GENCODE GTFs to not have this.

If your GTF does not contain gene lines, but can create them by merging all the exons in all the transcripts and filling in the gaps using `gtf2gtf` from `cgat`.

```
 cgat gtf2gtf --method=merge-transcripts --stdin Homo_sapiens.GRCh38.93.gtf.gz \
 | bedtools slop -b 10000 -g contigs.tsv -i - \
 | bedtools intersect -u -a intervals.bed -b -

```

# Answer
> 3 votes
First I prepared a bed file in which the gene intervals are augmented by 1KB before and after the gene start and end coordinates. Then I intersected this bed file with my original one with the option `-wa`, therefore retaining only the intervals in my original bed file that intersect with the bedfile produced from the original gtf, filtered by gene regions. The `-u` flag is also used to print only one line per overlap between the two files, in case an interval in the first overlaps with multiple intervals in the second.

```
# download information about chromosomes
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes ./

# Original GTF has a different notation than UCSC
cat Homo_sapiens.GRCh38.93.gtf | \
  sed -e $'s/^./chr&/' -e $'s/^chrMT/chrM/' | \
  awk '$3=="gene"' | \
  bedtools slop -b 1000 -g hg38.chrom.sizes -i - | \
  sed -e $'s/\tgene_id.*gene_name/\t/' -e $'s/gene_source.*//' | \
  tr -d '"; ' | \
  awk -v OFS=$'\t' '{print $1, $4, $5, $9}' > genes_1kb.bed

bedtools intersect -wa -u -a intervals.bed -b genes_1kb.bed > intervals.1kb.bed

```

---
Tags: gene, bed, gtf, interval
---


# Question
Title: Mapping RefSeq-ID or NCBI-ID (Entrez-ID) to EC number
I annotated my bacterial genomes using the new NCBI Prokaryotic Genome Annotation Pipeline and now, I want to annotate EC-numbers.

In the master annotation file (.gff), I get the RefSeq-ID of every successful annotation, e.g. `RefSeq:YP_805528.1`.

When I visit the respective NCBI page, I see that it sometimes is mapped to an EC number.

I'd like to automate this process: Convert every RefSeq-ID to an EC number. How can I do this?

(I may already have discovered the first step: I can map the RefSeq-ID to the NCBI-ID (Entrez-ID) using this NCBI conversion database: `YP_805528.1 -> 4419225`. But I don't see how to get from these databases to EC.)

Are there other ways? As far as I know, most people blast their genes to already-mapped databases like SwissProt and ExPASy ENZYME. But I don't like this approach because technically, the NCBI annotation pipeline already did this. Why compute the same thing twice?

# Answer
You can use Entrez Direct for this as follows:

```
efetch -db protein -id YP_805528.1 -format gpc -mode xml | xtract -insd Protein EC_number
YP_805528.1     6.3.2.13

```
> 2 votes

---
Tags: sequence-annotation, identifiers, phylogenetics, refseq
---


# Question
Title: levels of factors in the design have non-unique level names after make.names() is applied
When I try to import the data using `DESeq2` package:

```
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~Group + Cre)

```

I am getting the error:

> Error in validObject(.Object) : invalid class “DESeqDataSet” object: levels of factors in the design have non-unique level names after make.names() is applied. best to only uobject letters and numbers for levels of factors in the design

Here is how my `sampleTable` looks like:

```
    ID      Group Cre
1  Bla_36   IL4   -
2  Bla_38   IL4   -
3  Bla_11   IL4   -
4  Bla_17   IL4   -
5  Bla_20   IL4   +
6  Bla_21   IL4   +
7   Bla_8   IL4   +
8  Bla_22   IL4   +
9   Bla_7   IL4   +
10 Bla_35   PBS   -
11 Bla_37   PBS   -
12 Bla_13   PBS   -
13  Bla_1   PBS   -
14  Bla_9   PBS   +
15 Bla_23   PBS   +
16 Bla_19   PBS   +
17 Bla_24   PBS   +

```

What am I doing wrong here? And how to make it working?

# Answer
> 3 votes
There's a good reason to not allow `+` and `-` in the columns. It's because once you do that you get coefficients that can then start with a `+` or `-`, which is illegal in R and most other languages. In R, legal variable names consist of:

> A syntactically valid name consists of letters, numbers and the dot or underline characters and starts with a letter or the dot not followed by a number.

That doesn't really explain the error message you saw, of course. That can be solved by thinking about what R does to invalid character in names. Namely, it replaces them by a dot. So `+` and `-` get converted to `X.` and `X.`, which are of course identical and are what's really leading to the error message your saw. If you had used `+cre` and `-` it would have worked.

As an aside, you'll want to think in terms of positive and negative Cre if you talk to a biologist, since that's the appropriate terminology for it.

# Answer
> 2 votes
Weirdly enough, `+` and `-` are apparently not allowed in the column. After mapping them to longer strings it started to work:

```
library(plyr)
cre_id = revalue(samples$cre, c('-' = 'minus', '+' = 'plus'))

sampleTable <- data.frame(
  ID = samples$id,
  Group = samples$group,
  Cre = cre_id
)

```

and it looks like that:

```
       ID   Group   Cre
1  Bla_36   IL4 minus
2  Bla_38   IL4 minus
3  Bla_11   IL4 minus
4  Bla_17   IL4 minus
5  Bla_20   IL4  plus
6  Bla_21   IL4  plus
7   Bla_8   IL4  plus
8  Bla_22   IL4  plus
9   Bla_7   IL4  plus
10 Bla_35   PBS minus
11 Bla_37   PBS minus
12 Bla_13   PBS minus
13  Bla_1   PBS minus
14  Bla_9   PBS  plus
15 Bla_23   PBS  plus
16 Bla_19   PBS  plus
17 Bla_24   PBS  plus

```

I do not understand though, why...

---
Tags: r, deseq2
---


# Question
Title: Problem with SPAdes installation - no bin folder
I have been trying to download SPAdes for an important project I wanted to do, however, I have encountered a small problem while doing so. It appears that there is no 'bin' folder in my SPAdes installation file.

I have used 7-Zip as my main method of unpacking the file (As I am using Windows, I had to download the actual source code rather than a nice executable), and it came in the `.tar.gz` format, which I believe 7-Zip is rather good with.

A few pictures detailing my problems are shown below.

I hope I can get a little help on this matter. It is slightly time-sensitive, so I would really appreciate a little bit of assistance shortly. Thank you in advance!

As seen above, the first image is the command I used, the usual test command, the second is the file structure I got (Notice the absence of the 'bin' file), and the third is the error I received from running the `spades.py --test` command.

# Answer
> 1 votes
Simply unpacking the source code won't actually install much, you'll need to follow its installation instructions. Note that this might be possible on Windows if you use new bash command line, but in general don't expect any bioinformatics tools to work on windows.

---
Tags: assembly
---


title: List of human transcription factors

text:
# Question
Title: List of human transcription factors
I'm looking for a list of human transcription factor (TF) genes.

Where can I find a TF database or how can I retrieve this list from other databases?

# Answer
Lambert et al. (2018, Cell. 172(4):650-665) provide a manually curated list of human TFs at this website.
> 4 votes

---
Tags: database, data-retrieval, human
---


# Question
Title: Build a Pubmed query given long gene list
I would like to search the literature using the Pubmed search engine for a set of genes that are associated with a brain region of interest. Unfortunately, there are a lot of them. It's possible there is no solution if there is a character limit on the query. I think it would be easiest to do this with 'nix tools but I'm open to R libraries or Python implementations too.

The list is 908 genes long. It's from the Gensat annotation search engine, with the filters `Structure IS olfactory bulb` & `Expression level IS moderate to strong signal` & `Expression pattern IS region-specific`.

The list looks like, scraped from the dropdown list in the HTML source with a regular expression:

```
Aatf
Abcb1b
Abhd3
Abhd4
Abhd5
...

```

The Pubmed query should look like:

```
(((((((Aatf) OR Abcb1b) OR Abhd3) OR Abhd4) OR Abhd5) OR Abl2) AND olfactory bulb

```

but with 908 genes separated by OR and parentheses.

Bonus points if you find/write a utility which adds gene synonyms to the query!

# Answer
> 3 votes
After getting it out on paper (so to speak) I was able to accomplish what I wished with bash:

```
#!/bin/bash

gene_list=($(cat ./markers_clean.txt))
query=${gene_list[@]:0:1}
parens="("
for gene in ${gene_list[@]:1}
    do
        query="${query} OR ${gene})"
        parens="${parens}("
    done
echo "${parens}${query} AND olfactory bulb)"

```

With the output being:

```
((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((Aatf OR Abcb1b) OR Abhd3) OR Abhd4) OR Abhd5) OR Abl2) OR Ablim2) OR Acaa2) OR Ackr3) OR Acsbg1) OR Acvr1c) OR Acvr2a) OR Adam11) OR Adam23) OR Adcy1) OR Adcyap1) OR Adgrg1) OR Adora2a) OR Adra2a) OR Adra2c) OR AF529169) OR Agrp) OR Agt) OR Agtr1a) OR Agtr2) OR Ak4) OR Akap12) OR Akap5) OR Akr1c18) OR Aldh1l1) OR Aldh7a1) OR Aldoc) OR Als2) OR Amfr) OR Amigo1) OR Amigo2) OR Angpt2) OR Ank) OR Ank2) OR Ankrd34b) OR Api5) OR Apod) OR Arc) OR Arf6) OR Arhgap33) OR Arhgdig) OR Arhgef6) OR Arhgef7) OR Arl6ip1) OR Arpp21) OR Arvcf) OR Arx) OR Ascl1) OR Aspm) OR Atg4c) OR Atic) OR Atoh1) OR Atp1a1) OR Atp1b1) OR Atp2a2) OR Atp5a1) OR Atp5j2) OR Atp6v1b2) OR Atpif1) OR Avp) OR Avpr1a) OR Avpr1b) OR B2m) OR B3gat1) OR Bahcc1) OR Baiap2) OR Bcar1) OR Bdnf) OR Becn1) OR Bend5) OR Bmp1) OR Bmp6) OR Bmp7) OR Bnip3) OR Bsx) OR Btbd17) OR Btg1) OR Btg2) OR Btg3) OR C1ql2) OR C1ql3) OR C2cd2l) OR C33) OR Cacna1h) OR Cacna1i) OR Cacnb1) OR Cacng3) OR Cacng5) OR Cadps) OR Calb1) OR Calb2) OR Calca) OR Calcb) OR Calcr) OR Camk2a) OR Camkk1) OR Camkk2) OR Camkv) OR Card6) OR Cartpt) OR Casp9) OR Casr) OR Cbln1) OR Cbln2) OR Cbln4) OR Cbx5) OR Ccdc163) OR Ccdc184) OR Ccdc187) OR Ccdc3) OR Cck) OR Cckbr) OR Ccna2) OR Ccnd1) OR Ccne1) OR Ccne2) OR Ccsap) OR Cd63) OR Cdc42bpg) OR Cdc42ep4) OR Cdc6) OR Cdh1) OR Cdh11) OR Cdh3) OR Cdh6) OR Cdh9) OR Cdhr1) OR Cdk17) OR Cdk4) OR Cdk5) OR Cdk5r1) OR Cdk5r2) OR Cdk6) OR Cdkn1a) OR Cdkn1b) OR Cdkn2a) OR Cdkn2c) OR Cdkn2d) OR Cdon) OR Cdr2) OR Celf6) OR Cemip) OR Chd7) OR Chn2) OR Chrm1) OR Chrm4) OR Chrna2) OR Chrna3) OR Chrna4) OR Chrna5) OR Chrna7) OR Chrnb2) OR Chrnb4) OR Chst15) OR Chst2) OR Cic) OR Cisd1) OR Cited1) OR Ckmt1) OR Clcn1) OR Cldn1) OR Clu) OR Cnr1) OR Cntfr) OR Cntn2) OR Cntnap1) OR Coch) OR Col6a1) OR Colgalt2) OR Cpe) OR Cpne6) OR Crabp1) OR Creb1) OR Crh) OR Crhr1) OR Crhr2) OR Crim1) OR Crip2) OR Crk) OR Crym) OR Csnk1g3) OR Ctgf) OR Ctxn1) OR Cx3cl1) OR Cxcl11) OR Cxcl12) OR Cxcr4) OR Cxxc4) OR Cyb561) OR Cyp19a1) OR Cyp4v3) OR Cyr61) OR D8Ertd738e) OR Dbi) OR Dbp) OR Dcdc2a) OR Dct) OR Dctn1) OR Dcx) OR Ddit4l) OR Dgcr14) OR Dgcr6) OR Dhcr24) OR Dhrs3) OR Dhx57) OR Dkk3) OR Dlg4) OR Dlk1) OR Dll1) OR Dlx1) OR Dlx2) OR Dlx5) OR Dlx6) OR Dmbx1) OR Dnajb1) OR Dnajc12) OR Dnmt1) OR Dnmt3a) OR Doc2g) OR Drd1) OR Drd2) OR Drd3) OR Dsp) OR Dtx4) OR Dtymk) OR Dusp1) OR Dusp14) OR Dusp16) OR Dusp3) OR Dusp5) OR Dyrk1a) OR Dyx1c1) OR Dzip1) OR E2f8) OR Ebf2) OR Ebf3) OR Eef1a1) OR Efna1) OR Efnb1) OR Efnb2) OR Efr3a) OR Efs) OR Egr1) OR Eif1ax) OR Eif3h) OR Eif4e3) OR Eif4g2) OR Elavl2) OR Elavl3) OR Elavl4) OR Elfn1) OR Ell3) OR Emp1) OR Emx1) OR En2) OR Enah) OR Enc1) OR Eno1) OR Eno2) OR Entpd4) OR Eomes) OR Epha4) OR Epha7) OR Epha8) OR Eps8l2) OR Erp29) OR Esr2) OR Esyt3) OR Evx1) OR F2r) OR Fam171b) OR Fam81a) OR Fam84a) OR Fbxl5) OR Fbxo32) OR Fezf1) OR Fezf2) OR Fgd5) OR Fgf1) OR Fgf1) OR Fgf13) OR Fgf15) OR Fgfr1) OR Fign) OR Fkbp2) OR Fkbp3) OR Flrt1) OR Flrt2) OR Fn1) OR Fos) OR Foxi1) OR Foxi2) OR Foxm1) OR Foxp4) OR Frat1) OR Fryl) OR Frzb) OR Fxyd6) OR Fxyd7) OR Fzd2) OR Fzd3) OR Fzd5) OR Fzd8) OR Gaa) OR Gabbr1) OR Gabrb2) OR Gabre) OR Gabrg1) OR Gabrg2) OR Gad2) OR Gal) OR Galnt1) OR Galnt5) OR Galnt9) OR Galr3) OR Gap43) OR Gas1) OR Gbx2) OR Gck) OR Gcm2) OR Gdf1) OR Gfap) OR Gja1) OR Gjd2) OR Glce) OR Gldc) OR Glrb) OR Gls) OR Glud1) OR Gng1) OR Gng13) OR Gng3) OR Gng4) OR Gng5) OR Gng7) OR Gpc4) OR Gpm6a) OR Gpr12) OR Gpr139) OR Gpr151) OR Gpr165) OR Gpr21) OR Gpr27) OR Gpr37) OR Gpr68) OR Gpr83) OR Gpr88) OR Gprin3) OR Gpsm1) OR Gria2) OR Grid2ip) OR Grin2c) OR Grm1) OR Grm2) OR Grm4) OR Grp) OR Gsx2) OR Gtf2i) OR Gucy1a3) OR Hap1) OR Hbegf) OR Hcfc1) OR Hcn2) OR Hcn4) OR Hdac11) OR Hdc) OR Herc3) OR Hes6) OR Hey1) OR Hey2) OR Hgs) OR Hmgn2) OR Hmgn3) OR Hnrnpab) OR Hpcal1) OR Hrasls) OR Hrh2) OR Hsd11b1) OR Hspa12a) OR Hspb8) OR Htr1a) OR Htr1b) OR Htr1d) OR Htr3a) OR Htr3b) OR Htr4) OR Htr5b) OR Htr6) OR Icam5) OR Ict1) OR Id1) OR Id4) OR Idh3a) OR Ier5) OR Igf1) OR Igfbp2) OR Igfbp4) OR Igfbp5) OR Igsf9) OR Inpp5j) OR Insig2) OR Iqsec3) OR Irs4) OR Isl2) OR Islr2) OR Itgb3) OR Itm2b) OR Jag1) OR Jag2) OR Jak2) OR Jam3) OR Jph4) OR Jup) OR Kcnc1) OR Kcnc4) OR Kcnip2) OR Kcnj5) OR Kcnk1) OR Kcnq2) OR Kctd13) OR Kctd4) OR Kif11) OR Kif26a) OR Kit) OR Kitl) OR Klf1) OR Klf13) OR Kremen1) OR Krt19) OR Lamb1) OR Lamp5) OR Layn) OR Lbh) OR Ldb1) OR Lef1) OR Lfng) OR Lgals3) OR Lgi1) OR Lhx1) OR Lhx2) OR Lhx4) OR Lhx5) OR Lhx6) OR Lhx8) OR Limk1) OR Lix1) OR Lmo1) OR Lmo4) OR Lpar1) OR Lpl) OR Lrfn5) OR Lrig1) OR Lrig2) OR Lrp1) OR Lrp12) OR Lrp2) OR Lrrk2) OR Lrrn1) OR Lrrn3) OR Lrrtm2) OR Lss) OR Ly6h) OR Lynx1) OR Lypd1) OR Lypd2) OR Lypd6) OR Lzts1) OR Mafb) OR MAGEA8) OR Man1a2) OR Map1b) OR Map1lc3a) OR Map2k6) OR Map4k5) OR Mapk1) OR Mapk8) OR Mapt) OR March4) OR Matk) OR Mb21d2) OR Mbd2) OR Mbd4) OR Mc3r) OR Mc4r) OR Mdk) OR Med12) OR Mef2c) OR Meis1) OR Met) OR Mfap2) OR Mfap4) OR Mfn1) OR Mmp17) OR Mmp2) OR Mpc1) OR Mst1r) OR Msx1) OR Mtap) OR Mturn) OR Mxi1) OR Mycn) OR Nae1) OR Napb) OR Nars) OR Ncor1) OR Ndn) OR Ndnf) OR Ndrg1) OR Ndufv1) OR Ndufv2) OR Nectin3) OR Nefh) OR Nefl) OR Nefm) OR Neurod1) OR Neurod4) OR Neurod6) OR Neurog1) OR Neurog2) OR Nfil3) OR Nhlh1) OR Nhp2) OR NKAIN2) OR Nkain3) OR Nkx2) OR Nlgn3) OR Nmb) OR Nmbr) OR Nmu) OR Nnat) OR Nog) OR Nomo1) OR Nos1) OR Notch1) OR Notch2) OR Nov) OR Nova1) OR Npc1) OR Npc2) OR Nptx1) OR Nptx2) OR Npy) OR Npy1r) OR Npy2r) OR Npy5r) OR Nr2f2) OR Nr4a1) OR Nrep) OR Nrgn) OR Nrip3) OR Nrp1) OR Nrtn) OR Nrxn2) OR Nsf) OR Nsg1) OR Nsmf) OR Ntn1) OR Ntng2) OR Ntrk1) OR Ntrk2) OR Ntrk3) OR Nts) OR Nudt2) OR Numb) OR Numbl) OR Odc1) OR Olfm1) OR Olfr734) OR Olig2) OR Oma1) OR Ophn1) OR Opn3) OR Opn4) OR Oprk1) OR Otof) OR Otx2) OR Oxt) OR Oxtr) OR P2rx7) OR Pacsin3) OR Pak1) OR Pard6a) OR Pcbd1) OR Pcdh2) OR Pcdh9) OR Pcnt) OR Pcp4) OR Pcsk9) OR Pde1b) OR Pde1c) OR Pdp1) OR Pdyn) OR Peg3) OR Penk) OR Per2) OR Pex5l) OR Pgrmc1) OR Phactr1) OR Phactr4) OR Phox2b) OR Pias3) OR Pik3r1) OR Pirt) OR Pkd1) OR Plagl1) OR Plcxd2) OR Pld3) OR Plekhg1) OR Plk2) OR Plk3) OR Plp1) OR Plxdc1) OR Plxna1) OR Plxnc1) OR Pmp22) OR Pnmt) OR Pnoc) OR Podxl2) OR Polr1d) OR Polr3e) OR Pomc) OR Pou2f1) OR Pou3f4) OR Ppfia4) OR Ppm1e) OR Ppp1ca) OR Ppp1r14c) OR Ppp1r16b) OR Ppp1r17) OR Ppp1r1a) OR Ppp1r1b) OR Ppp1r2) OR Ppp1r7) OR Ppp2r2a) OR Ppp4c) OR Ppp6c) OR Prkacb) OR Prkcd) OR Prkcq) OR Prmt2) OR Prnp) OR Prodh) OR Prok2) OR Prokr2) OR Prox1) OR Prss12) OR Prtg) OR Psen1) OR Psen2) OR Ptch2) OR Pten) OR Ptf1a) OR Ptgds) OR Ptn) OR Ptpn2) OR Ptpn5) OR Ptprf) OR Pus1) OR Pygb) OR Rab37) OR Rac1) OR Ramp3) OR Rap1a) OR Rasa3) OR Rasal1) OR Rasgrp1) OR Rasgrp2) OR Rasl11b) OR Rbp1) OR Rcan2) OR Rcn1) OR Rel) OR Rest) OR Rfng) OR Rgs12) OR Rgs16) OR Rgs2) OR Rgs3) OR Rgs4) OR Rgs5) OR Rgs8) OR Rgs9) OR Rhob) OR Rhou) OR Rif1) OR Rlbp1) OR Rnd1) OR Rnd2) OR Rnd3) OR Rorb) OR Rpl18) OR Rpl31) OR Rprm) OR Rraga) OR Rtn4rl1) OR Rubcn) OR Rufy3) OR Runx1) OR Rxfp3) OR S1) OR Satb2) OR Scg2) OR Scg3) OR Scn3b) OR Scrn1) OR Scrt1) OR Scube1) OR Scube2) OR Sdc3) OR Selenom) OR Selenow) OR Sema3a) OR Sema3b) OR Sema3c) OR Sema3e) OR Sema3f) OR Sema4a) OR Sema6d) OR Sema7a) OR Sept3) OR Sept6) OR Serf1) OR Serpina9) OR Serpini1) OR Sertm1) OR Sez6) OR Sez6l2) OR Sf1) OR Sfrp1) OR Sfrp2) OR Sfrp5) OR Sim1) OR Sirpa) OR Sirt2) OR Skil) OR Sla) OR Slc1) OR Slc12a5) OR Slc13a3) OR Slc17a6) OR Slc17a7) OR Slc18a2) OR Slc1a1) OR Slc1a2) OR Slc1a3) OR Slc1a6) OR Slc22a3) OR Slc25a12) OR Slc32a1) OR Slc33a1) OR Slc41a3) OR Slc6a1) OR Slc6a12) OR Slc6a2) OR Slc6a3) OR Slc6a6) OR Slc6a9) OR Slc7a14) OR Slc8a2) OR Slc9a3r1) OR Slit1) OR Slitrk2) OR Slitrk3) OR Slitrk4) OR Slitrk5) OR Slitrk6) OR Smad6) OR Smarce1) OR Smoc2) OR Snap25) OR Snap91) OR Sncg) OR Snhg11) OR Snx17) OR Sort1) OR Sox1) OR Sox11) OR Sox2) OR Sox4) OR Sox8) OR Sox9) OR Sp4) OR Sprn) OR Spry1) OR Sptan1) OR Sqle) OR Srgap3) OR Srsf7) OR Sstr1) OR Sstr2) OR St8sia5) OR Stard5) OR Stau1) OR Strip2) OR Strn) OR Stum) OR Stx1a) OR Sufu) OR Sulf2) OR Sun2) OR Suv39h1) OR Sv2a) OR Sv2b) OR Syndig1l) OR Syngr3) OR Synj2) OR Syt13) OR Syt17) OR Syt2) OR Syt3) OR Syt4) OR Syt6) OR Syt7) OR Tac1) OR Tac2) OR Tacr3) OR Taf13) OR Tagln3) OR Tbc1d8) OR Tcf3) OR Tgif1) OR Thbs1) OR Thbs2) OR Thbs3) OR Thbs4) OR Them6) OR Thra) OR Thrsp) OR Tigar) OR Tle1) OR Tle3) OR Tmeff1) OR Tmem35a) OR Tmem59l) OR Tnc) OR Tor1b) OR Tpbg) OR Tpd52l1) OR Tph2) OR Tppp3) OR Traf4) OR Trappc3) OR Trf) OR Trh) OR Trhr) OR Trib2) OR Tsc22d2) OR Tspan14) OR Tspan7) OR Tuba1b) OR Tubb3) OR Tubb4a) OR Ube3a) OR Ubxn11) OR Uchl1) OR Ucp2) OR Uhmk1) OR Uncx) OR Uso1) OR Vat1) OR Vcan) OR Vdac2) OR Vegfa) OR Vim) OR Vip) OR Vipr2) OR Vldlr) OR Vmn1r2) OR Vsnl1) OR Vsx2) OR Washc2) OR Wdr46) OR Wfs1) OR Wif1) OR Wnt2b) OR Wnt5b) OR Wnt7a) OR Wnt8b) OR Wrb) OR Wscd1) OR Wwc1) OR Ypel4) OR Zbtb18) OR Zbtb2) OR Zcchc12) OR Zfp5) OR Zfp532) OR Zic1) OR Zic2) OR Zic3) OR Zic4) OR Zic5) OR Zwint) AND olfactory bulb)

```

As far as I can tell, this does not surpass any character limits on Pubmed.

---
Tags: text-processing, literature-search
---


title: Codon usage analysis for whole genomes

text:
# Question
Title: Codon usage analysis for whole genomes
I am new to bioinformatics. So if these questions seem you to a bit childish please forgive me.

I have two queries.

1. I am intending to perform a codon usage analysis followed by correspondence analyses for multiple microbial whole genomes of one bacterial species to find the association with the isolation source and ST type. Is this the right approach?
2. I have concatenated all the CDS in a single genome and joined all the genomes (average genes per genome is 4100). The problem is all the CDS when concatenated also have the stop codons in it. CodonW, famous for this analyses, cannot begin the analyses in the presence of stop codons. How can I remove stop codons?

Kindly suggest me a solution to the problem and concept.

# Answer
This python script will remove all stop codons from your fasta file. I called it `remove_stops.py`

```
#!/usr/bin/env python3
import sys
from Bio import SeqIO
stop_codons = ["TAG", "TGA", "TAA", "UAG", "UGA", "UAA"]

fasta_it = SeqIO.parse(open(sys.argv[1]), 'fasta')
for fasta in fasta_it:
    name, sequence = fasta.description, str(fasta.seq)
    try:
        assert(len(sequence) % 3 == 0)
    except:
        raise IOError("length of sequence {} was not divisible by three".format(name))
    print(">{}".format(name))
    print("".join([sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3)
                   if sequence[i:i+3] not in stop_codons]))

```

You can run it like so:

```
python remove_stops.py stop_codons.fasta > filtered.fasta

```

This program won't work if there are any fasta records that do not have a length divisible by three. This indicates a sliding translation stop and these sequences should be trimmed individually. Also, all the genes must be in the 5'-\>3' direction.

It works by looking at each codon `[sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3)`, then just makes sure they are not a stop codon: `if sequence[i:i+3] not in stop_codons]`
> 2 votes

---
Tags: sequence-analysis, software-usage, codon
---


title: Replace lowercase characters with -

text:
# Question
Title: Replace lowercase characters with -
I have an output from vcfutils.pl vcf2fq with specified minimal depth, and it means that nucleotides with not enough depth are lowercase. I would like to change them to gaps. I could do it in higher language but I would like to know the solution in bash. I have found simple solution to replace all lowercase to uppercase:

```
awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' in.fna > out.fna

```

but I would like to replace lowercase with -. The first row need to be omitted. Could you help me? Thanks

# Answer
The following `sed` command will do the trick.

`sed -e '/^[^>]/ s/[a-z]/\-/g' in.fna > out.fna`

Lines starting with `>` will be escaped and the lowercase characters will be replaced with `-` .
> 4 votes

# Answer
You could try it with `sed`.
Replace all lower-case chars with `-` in lines not starting with `>`:

```
sed -e '/^>/!s/[a-z]/-/g' in.fa

```
> 3 votes

# Answer
You can also use this solution where **tr** command from linux can solve your problem:

`tail -n +2 in.fna | tr '[:lower:]' '-'`
> 1 votes

---
Tags: bash, shell
---


# Question
Title: How do I interpret my phylogenetic tree construction from Glucosidases in the Actinobacteria class?
Dear knowledgeable bioinformaticians, I'm currently pursuing an undergraduate degree and as a result I do not have so much knowledge in this topic. I'm having some troubles understanding my phylogenetic three construction.

For the construction of a phylogenetic tree UGENE was used. 11 homologous proteins from different species in the same genus/phyla were aligned. For the tree construction PHYLYP Neighbour Joining (NJ) with the distance matrix model “Jones-Taylor-Thornton” was used. Two trees were constructed, one without bootstrapping and one with bootstrapping (1000 pseudoreplicates).

How do I interpret my values? Are they unusable since I have low bootstrap-values in a lot of my nodes? All guidance is appreciated, thank you.

# Answer
Bootstrap values usually give you useful information, be they low or high. Usually, a "high" bootstrap value is something \>= 70%. That means a clade is strongly supported; lower values are telling you the clade is not clearly separated from others. This could be for several reasons, e.g.:

1) A subpopulation is not really separated from others, e.g. the viral sequences in town "A" and town "B" are completely admixed, meaning the population is panmictic and there is continuous gene flow between the two of them;

2) the populations diverged not long ago, it's too soon to have any phylogenetic signal, at least using the marker you are using (different genes evolve at different rates). Usually, housekeeping genes (which are usually very conserved) are used for populations that supposedly diverged long time ago.

The interpretation of your tree should take into account the biology and epidemiology of your target organism. If you really want to learn more, I suggest a lot of reading, and definitely a class dedicated to molecular phylogenetics.

This might be a good start: https://www.amazon.com/Phylogenetic-Handbook-Practical-Approach-Hypothesis/dp/0521730716

Also, the NJ algorithm should allow for uneven branch lengths, maybe the tree you are showing is the consensus tree?

Hope it helps

Max
> 2 votes

---
Tags: phylogeny, figure-explanation, phylogenetics
---


# Question
Title: Using ONT MinION, why is that we cannot get a full length DNA read?
Here, https://www.youtube.com/watch?v=E9-Rm5AoZGw, ONT claims they can read non-fragmented DNA and that with the hairpin like structure the double strands can be read at once. Then why is that MinION reads are of variable in length and not the length of a whole DNA? BUt here:https://www.youtube.com/watch?v=e9wtAIvPPxY, ONT claims that the read lenght can be configured, what does this mean?

# Answer
> 2 votes
The MinION should be able to sequence DNA fragments entirely, regardless of the size, if you can get them to the pore intact. That is the problem: during DNA extraction and library preparation you will introduce breakage.

---
Tags: nanopore, long-reads, minion
---


# Question
Title: What does that mean if there is no bifurcations in my time series single cell data?
Sorry,

I have analysed my time series single cell data with Scuba algorithm but that says that there is no bifurcations in my data and outputted a single chain of 4 states with no branching like below figure

Sorry, I googled but I am not quite sure what this means

Any suggestion please?

Thanks a lot

# Answer
> 7 votes
The most obvious interpretation is that the four cell lines are related in a linear pathway rather than a branching one.

---
Tags: single-cell, matlab
---


# Question
Title: Collapse cell barcodes distribution within 1 Hamming distance
I have a barcode distribution from single-cell data, e.g:

```
11612552 TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
...

```

where the first column is the number of reads associated with the cell barcode in the second column.

I would like to collapse each set of cell barcodes within 1 Hamming distance (i.e. allowing one mismatch between barcodes), keeping only the barcode with the highest read count and adding the read counts of the duplicates to it.

I was thinking to use `awk` and `grep` to filter the barcode list in an iterative procedure using regular expressions, but I think this would require too much time. I can share my code if needed.

Is there a way to do it using unix tools or already published pipelines (e.g. UMI tools)?

This is my code at the moment (it is not optimized though):

```
#!/bin/awk -f

{
    matching = 0;
    for (r in regexs) {
        if ($2 ~ regexs[r]) {
        matching = 1;
        counts[r] += int($1);
        break;
    }
}
if (!matching) {
    regex = "";
    for (i = 1; i <= length($2); i++) {
        l_regex = substr($2, 1, i-1) "." substr($2, i+1);
        if (length(regex) == 0)
            regex = l_regex;
        else
            regex = regex "|" l_regex;
    }
    regexs[$2] = regex;
 counts[$2] = int($1);
    }
    # so we know that 100 records have been processed
    if (!(NR % 100)) printf "." > "/dev/stderr";
}
END {
    for (c in counts) {
        print counts[c], c | "sort -k1 -rn";
    }
}

```

# Answer
This doesn't use standard unix tools, but rather the API for UMI tools:

```
#!/usr/bin/env python
from umi_tools import network

cb = []  # cell barcodes
cbc = dict()  # cell barcode counts
for line in open("UMIs.txt"):
    cols = line.strip().split("\t")
    cols[1] = bytes(cols[1], "ascii")
    cb.append(cols[1])
    cbc[cols[1]] = int(cols[0])

uc = network.UMIClusterer()
CBclusters = uc(cb, cbc, 1)
cbFinal = dict()
for l in CBclusters:  # This is a list with the first element the dominant barcode
    cbFinal[l[0]] = 0
    for x in l:  # Iterate over all barcodes in a cluster
        cbFinal[l[0]] += cbc[x]
print(cbFinal)

```

If you give it input in `UMIs.txt` of the following form (note the last two lines differ from the first with a distance of 1):

```
11612552    TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
20  TCCAGAGCACTGCATAACTCAA
2   TCCTGAGCACTGCATAACTCNA

```

Then it will produce:

```
{b'TAAGGCGATCTCTCCGTACGTG': 5985219,
 b'TAAGGCGAACTGCATAACATGT': 6626630,
 b'GCTACGCTACTGCATAAGTCCA': 9349711,
 b'CGTACTAGTCTCTCCGCGGCTA': 8161950,
 b'TCCTGAGCGTAAGGAGCAGATC': 8102383,
 b'CAGAGAGGCTAAGCCTGCACAT': 8343678,
 b'TCCTGAGCACTGCATAACTCAA': 11612574,
 b'AGGCAGAATCTCTCCGAGTGCT': 6210446,
 b'AAGAGGCAGTAAGGAGATAATC': 7110298,
 b'CGAGGCTGTATCCTCTACTAGA': 6390489}

```
> 3 votes

# Answer
The difficulty with all these things comes with what to do about barcodes that are 1 edit distance away from more than one other barcodes in ways that form complex networks.

Cell barcodes are a bit different from UMIs. Partly this is because we expect them to be sparser. Partly because the consequences of inncorect error correction at the CB level are much worst than at the UMI level. The way that UMI-tools does its Cell Barcode error correction is:

* Determine a whitelist of barcodes using the "knee" method
* Iterate through the list of non-whitelist barcodes and check if they are within one edit distance of one **and only one** whitelist barcode. (this is done with `umi_tools.umi_methods.getErrorCorrectMapping` if you are interested).
* Things that pass the above are collapsed into the whilelist barcode, else they are discarded.

The function `umi_tools.umi_methods.getCellWhitelist` takes a list of barcodes and their frequencies and returns a whitelist and (optionally) a mapping of which barcodes should be corrected to which others.
> 3 votes

---
Tags: single-cell, barcode
---


# Question
Title: Retrieve all genes for a specific gene ontology
I'm having a bit of confusion about how can I retrieve all the genes that are associated to a specific gene ontology.

Lets say I need to find all the genes that are associated with this GO entry:

> GO:0030098 lymphocyte differentiation

This is the graph representing the parents of my GO. In the same page I find the **Child Terms** that represent the children of my term.

Now, I can download the goa file (related to the organism I'm interested into) containing the GO associations and `grep` my term to extract a list of associated genes but, as far as I can understand, this list does not represent the whole graph that actually describe my GO term.

Should I also get the genes in the childred terms of my term? Or in the parents terms? If yes, then again I have to go back in the graph and find the parents/children of the found parents/children and so on, till the whole graph is visited.

Is this the correct approach? Or should I only focus on those genes related to my GO:0030098, and that's it?

I found this: Retrieving a list of human genes having GO associations and this but I think they don't address my question since they do not describe the parent-children association.

Also I found this, on the GO FAQ, but what I need to do is quite the opposite.

Summarising: I need to find the genes associated to a specific GO term. Should I also look for its parents or children or not?

Any clues?

# Answer
You can use BioMart to filter by your GO term and get the genes as attributes. BioMart is ontology-aware so will pull out all genes associated with your term and with any of its child terms.

There is no need to look up the GO child terms as BioMart already deals with this. A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.

If you're not already familiar with BioMart, there's a help video.

Here's a BioMart query that gets you exactly what you need in human using the web interface

If you want to work with the R interface to biomaRt, here's some documentation and this is the identical query to above:

```
library(biomaRt)
ensembl=useEnsembl(biomart="ensembl", dataset='hsapiens_gene_ensembl')
getBM(attributes = c('ensembl_gene_id','external_gene_name', 'go_parent_term', 'name_1006'),
      filters = 'go',
      values = 'GO:0030098',
      mart = ensembl)

```
> 4 votes

---
Tags: gene-ontology
---


# Question
Title: Group genes by functional categories suming expression values
Using the count of rpkm values from genes in a metagenome sample, I want to group these genes into established categories (for example KEGG or COG). For each sample, my goal is to determine which categories are better represented in each one.

Considering the example table whit genes, samples (1-5) and categories. Each gene has their expression and the column categories classify them.

```
import sys
import numpy

dictionary = {}

table = open(sys.argv[1], "r")

for line in table:
    next()
    cols = line.strip().split("\t")
    gene = cols[0]
    s1 = float(cols[1])
    s2 = float(cols[2])
    s3 = float(cols[3])
    s4 = float(cols[4])
    s5 = float(cols[5])
    counts = [s1, s2, s3, s4, s5]
    kegg = cols[6].strip().split(" ")[0].replace('"','')

    if kegg not in dictionary:
        dictionary[kegg] = [counts]
    else:
        dictionary[kegg].append(counts)

#print dictionary
for k, v in sorted(dictionary.iteritems()):
    m = numpy.array(v)
    print k,  [sum(m[:,i]) for i in range(5)] # range(number_of_samples)

```

Original table to count the columns values:

```
gene    s1  s2  s3  s4  s5  Category
name01  0   2   2   0   0   A
name02  3   0   1   0   0   A
name03  0   2   1   0   0   A
name04  0   0   1   0   0   B
name05  5   0   1   0   0   C
name06  1   0   0   0   0   D
name07  2   0   0   0   0   D
name08  0   0   3   0   0   E
name09  1   0   0   0   0   F
name10  1   0   0   0   0   F
name11  3   0   0   0   0   F

```

The result obtained is this:

```
"""
# eliminated from the question

Cat s1  s2  s3  s4  s5
A   0.0 2.0 1.0 0.0 0.0
B   0.0 0.0 1.0 0.0 0.0
C   5.0 0.0 1.0 0.0 0.0
D   2.0 0.0 0.0 0.0 0.0
E   0.0 0.0 3.0 0.0 0.0
F   3.0 0.0 0.0 0.0 0.0
"""

```

But my expected (obtained) table is this:

```
Cat s1  s2  s3  s4  s5
A   3.0 4.0 4.0 0.0 0.0
B   0.0 0.0 1.0 0.0 0.0
C   5.0 0.0 1.0 0.0 0.0
D   3.0 0.0 0.0 0.0 0.0
E   0.0 0.0 3.0 0.0 0.0
F   5.0 0.0 0.0 0.0 0.0

```

# Answer
> 4 votes
As suggested in the comments, the pandas module is quite convenient for this type of work.

Here is how you could do (assuming you have tab-separated values as input, and you want similar file format for the output):

```
import pandas as pd

counts = pd.read_table("counts.tsv")
# At this point, counts looks as follows:
#       gene  s1  s2  s3  s4  s5 Category
# 0   name01   0   2   2   0   0        A
# 1   name02   3   0   1   0   0        A
# 2   name03   0   2   1   0   0        A
# 3   name04   0   0   1   0   0        B
# 4   name05   5   0   1   0   0        C
# 5   name06   1   0   0   0   0        D
# 6   name07   2   0   0   0   0        D
# 7   name08   0   0   3   0   0        E
# 8   name09   1   0   0   0   0        F
# 9   name10   1   0   0   0   0        F
# 10  name11   3   0   0   0   0        F

summed_by_cat = counts.groupby("Category").sum()
# At this point, summed_by_cats looks as follows:
#           s1  s2  s3  s4  s5
# Category
# A          3   4   4   0   0
# B          0   0   1   0   0
# C          5   0   1   0   0
# D          3   0   0   0   0
# E          0   0   3   0   0
# F          5   0   0   0   0

# Change the index name to get the desired first column header in the output:
summed_by_cat.index.name = "Cat"

# Write to a file
summed_by_cat.to_csv("by_cat.tsv", sep="\t")

```

And resulting the file contains the following:

```
Cat s1  s2  s3  s4  s5
A   3   4   4   0   0
B   0   0   1   0   0
C   5   0   1   0   0
D   3   0   0   0   0
E   0   0   3   0   0
F   5   0   0   0   0

```

You could also load your table using `counts = pd.read_table("data.tsv", index_col="gene")` and `counts` would have the gene names as index instead of default numbers:

```
        s1  s2  s3  s4  s5 Category
gene
name01   0   2   2   0   0        A
name02   3   0   1   0   0        A
name03   0   2   1   0   0        A
name04   0   0   1   0   0        B
name05   5   0   1   0   0        C
name06   1   0   0   0   0        D
name07   2   0   0   0   0        D
name08   0   0   3   0   0        E
name09   1   0   0   0   0        F
name10   1   0   0   0   0        F
name11   3   0   0   0   0        F

```

This doesn't seem to affect the final outcome.

# Answer
> 1 votes
Below is code to do that, though if these are really gene counts then I highly suggest you provide more details about why you're trying to do this before you actually do it.

```
#!/usr/bin/env python
print("Cat s1  s2  s3  s4  s5")
v = []
category = None
for line in open("counts.txt"):
    if line.startswith("gene"):
        continue
    cols = line.split()
    if cols[-1] == category:
        for idx in range(5):
            v[idx] += int(cols[idx + 1])
    else:
        if category:
            print("{}   {}".format(category, " ".join(["{:.1f}".format(x) for x in v])))
        v = [int(x) for x in cols[1:-1]]
        category = cols[-1]
print("{}   {}".format(category, " ".join(["{:.1f}".format(x) for x in v])))

```

---
Tags: python
---


# Question
Title: How to convert data in gmt format to dataframe?
I downloaded c5: gene ontology gene sets file from http://software.broadinstitute.org/gsea/downloads.jsp

I opened the "c5.all.v6.2.symbols.gmt" file in csv format and It looks like below:

I want to convert the .gmt file into a dataframe. And it should look like below:

can anyone say how I can do that. Thanq

# Answer
Ok. Answering my own question.

I did this way. And I got the output I want.

```
install.packages("msigdbr")
library(msigdbr)

m_df = msigdbr(species = "Homo sapiens", category = "C5")
head(m_df)

```

Output:

```
gs_name gs_id gs_cat gs_subcat human_gene_symb~ species_name entrez_gene
  <chr>   <chr> <chr>  <chr>     <chr>            <chr>              <int>
1 GO_14_~ M184~ C5     MF        AANAT            Homo sapiens          15
2 GO_14_~ M184~ C5     MF        AKT1             Homo sapiens         207
3 GO_14_~ M184~ C5     MF        ARRB2            Homo sapiens         409
4 GO_14_~ M184~ C5     MF        BAD              Homo sapiens         572
5 GO_14_~ M184~ C5     MF        DAB2IP           Homo sapiens      153090
6 GO_14_~ M184~ C5     MF        DDIT4            Homo sapiens       54541
# ... with 2 more variables: gene_symbol <chr>, sources <chr>

```

found this solution from here https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html
> 3 votes

# Answer
You can import the csv file as a table. This does not require further libraries:

```
#you're reading a csv file, using tab as field separator and considering the first line as the headers of the data table
data <- read.csv(filename, sep="\t", header=TRUE)

```

then what you need are just two columns of your **data**:

```
colnames(data) #will give you the names of the columns that can be accessed
#to access two columns and merge them into a selectedColumns table
selectedColumns <- cbind(data$colName1, data$colName2)

```
> 1 votes

# Answer
You can use GSA.read.gmt function from GSA package. The following code can be used to convert the file to a dataframe. Just ignore the warnings.

**Original\_response**

```
library(GSA)
data <-  GSA.read.gmt("c5.all.v6.2.symbols.gmt")
gene_names <- unlist(data$genesets, use.names=FALSE)
your_dataframe <- cbind(data$geneset.names,gene_names)
colnames(your_dataframe) <- c("Pathways","Genes")
mydataframe <- as.data.frame(your_dataframe)
head(mydataframe)
                                               Pathway   Gene
1            GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2C
2                           GO_CARDIAC_CHAMBER_DEVELOPMENT POLR2J
3 GO_DNA_DEPENDENT_DNA_REPLICATION_MAINTENANCE_OF_FIDELITY  CTDP1
4                                      GO_CIRCADIAN_RHYTHM   RDBP
5              GO_PHOSPHATIDYLSERINE_ACYL_CHAIN_REMODELING COBRA1
6                               GO_SPINAL_CORD_DEVELOPMENT   RSF1

```

**Edited\_response**

You can use the code given below to achieve your desired output:

```
library(GSA)
data <-  GSA.read.gmt("c5.all.v6.2.symbols.gmt")
len_vec=c()           # Now create a vector for containing the length of genes at each position
len_vec[1] = 3
for(i in 1:length(data$genesets)){len_vec[i] <- c(length(data$genesets[[i]]))}
pathway_vec <- unlist(Vectorize(rep.int)(data$geneset.names, len_vec),use.names = FALSE) # Now create a vector for all the pathways in the data
desired_df <- as.data.frame(cbind(pathway_vec,unlist(data$genesets,use.names = FALSE))) # This gives your desired dataframe
head(desired_df)
                                    Pathway      Gene
1 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2C
2 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION POLR2J
3 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION  CTDP1
4 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION   RDBP
5 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION COBRA1
6 GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION   RSF1

```
> 1 votes

---
Tags: r, file-formats, format-conversion, gsea, conversion
---


title: Error correction within the long read

text:
# Question
Title: Error correction within the long read
I am going to get some data from plasmid sequencing to identify SNPs on the plasmids. What is done in the lab is the following:

* The plasmids are purified by size.
* We amplify the plasmids using the phi29 polymerase. The polymerase will go through the plasmid multiple times. Hence we get the same sequence concatenated multiple times. **My question is related to this step:**
* Finally, We sequence it using Oxford Nanopore.

My question is:

On step two, I wrote that "we get the same sequence concatenated multiple times". For me, this a potential source of information to correct the base calls prior to aligning them to the reference. Since you have the same sequence multiple times (concatemers).

What I would like to know is if there is a tool that uses the information of the concatemers to improve the base calls. I have tried to find some methods but found none. I know that PacBio has a similar flavour using the "circular consensus calling", but I have not found any methodological explanation.

thanks!

# Answer
Yes, there is a tool to do this called R2C2 by the Vollmers lab at UC Santa Cruz.

https://www.biorxiv.org/content/early/2018/06/04/338020
> 2 votes

---
Tags: nanopore, long-reads, minion, literature-search, minimap2
---


# Question
Title: Accessing UCSC genome via ssh results in a validation error
I am working on macOS High Sierra.

I am following the steps described here to lift an annotation over from one version of a genome to another. I am now just using the example genomes provided in the tutorial.

At this step:

```
time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

```

I get the following error in the `do.log` file:

```
Warning: Permanently added 'localhost' (ECDSA) to the list of known hosts.
username@localhost: Permission denied (publickey,password,keyboard-interactive).
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Thu Nov  8 14:27:32 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
align: localhost does not have /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.ooc -- if that is not the correct location, please run again with -ooc.

```

I think the main problem is here:

```
username@localhost: Permission denied (publickey,password,keyboard-interactive).

```

This post describes this issue. However:

* I do not have this file: ~/.ssh/authorized\_keys. My files in ~/.ssh are: `github_rsa github_rsa.pub known_hosts`. One of the lines of `known_hosts` is localhost. Does this presumably allow me to ssh to my own laptop? I still don't understand why I need to ssh in this script if it's run locally.
* `PasswordAuthentication yes` is already in `/etc/ssh/sshd_config`
* If I run the command without the `time` command, I still get the same error. I am not prompted to give a password.

Could anyone help me overcome this problem, please?

# Answer
> 2 votes
You need to setup ssh keys. All the commands in the script are going to use ssh to run the commands. Note the tips here of how to setup ssh keys:

http://genomewiki.ucsc.edu/index.php/Parasol\_job\_control\_system#SSH\_keysssh key setup

---
Tags: sequence-annotation, ucsc, errors, liftover
---


title: Finding URL for UCSC data hub

text:
# Question
Title: Finding URL for UCSC data hub
I am trying to build a UCSC track hub (a private one).

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

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

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

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

```

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

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

```

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

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

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

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

```
<html><head><title>my home directory</title></head>
<body><h1>hello world</h1></body></html>

```

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

---
Tags: ucsc
---