title: Reading multiple raw files in Seurat

text:
# Question
Title: Reading multiple raw files in Seurat
I have multiple single cell samples to analyze and I'm following the instructions in Satija Lab's website. I want to merge all the count files from all the samples at once, and associate the metadata to each sample.

For the `Read10x` command, is there a way to read multiple files at once, such as merging all the filter\_featire\_bc\_matrix files?

I have all the samples saved in a directory on my desktop, is it possible to call the directory within the script for it to read all the raw counts?

This is what I have thus far:

```
# load the required libraries
library(dplyr)
library(Seurat)
library(ggplot2)

-----------------------------------------------

# set up the working directory

wd = ("~/Seurat") setwd(wd)

# read all metadata

metadata = read.table("metadata_Jan2020.csv", sep=';', header=T, stringsAsFactors=F)
rownames(metadata) = metadata$LibraryID

#metadata
#mf = metadata

mf = metadata %>% filter(Type %in% c("PDAC", "NORMAL"))
mf = mf %>% filter(!(ID %in% c("PDAC-HTB2854")))`

```

# Answer
Seurat has a vignette specifically for combining multiple 10x libraries.
> 2 votes

---
Tags: r, scrnaseq, seurat, 10x-genomics
---


title: What does this statement mean?

text:
# Question
Title: What does this statement mean?
I have a copy number segment file with this by SCAT R package

```
Chromosome  Start   End
1           10583   120983
1           121009  245455
1           404187  5719467
1           5719471 5735807

```

I want to count the number of markers in these locations

I asked somebody how to do this, he replied:

**We created marker file by dividing the genome into 10 kb bins, so the markers file contains the start and stop positions for each of these 10kb segments for each chromosome along with ASCAT segments start and end positions. Then, for each ASCAT copy number segment, we count how many markers fall within that region for total probes.**

Can somebody helps to understand what he mean by dividing the genome

# Answer
> 4 votes
They created 10kb bins, which could be called "dividing the genome". These are just 10kb long contiguous stretches (so position 1-10000,then 10001-20000, etc.). This is, of course, a nonsensical way of going about things, because "counting those up" in the regions you posted is the same as taking the width of the regions, dividing by 10000, and rounding a bit.

Presumably you have your markers in something like BED format, so just use bedtools to intersect them and do the counting for you.

---
Tags: r, wgs, cnv
---


# Question
Title: What are the different kinds of bioluminescent genes?
I know of the common green glow gene but I forgot the name and I also know that some algae glow blue. There are so many types of bioluminesent organisms, so I am wondering what species have which genes and what other genes are associated with their metabolism and the systems that regulate these genes.

# Answer
> 6 votes
There are three processes wherein light is emitted.

* **Bioluminescence**: a chemical reaction releases light. The enzyme that does this is luciferase, while its substrate is luciferin, a small molecule. There are half a dozen organisms that have analogous (evolved independently) enzymes and totally different substrates, such as fireflies, Alivibrio bacteria, diatomes, deep sea fish, some mushrooms etc. The wiki article has a list of these and their small molecules.
* **Fluorescence**: light is absorbed by the molecule and is emitted at a lower energy (wider wavelength, redshifted). Green fluorescent protein of Aequorea victoria jellyfish is the classic example of this. It is excited with UV and emits green light. Here is an interactive model, I threw together. Variants exists for biotech reasons that span the whole of the visible spectrum. Other species have a homologous beta barrel, such as corals (red).
* **Phosphorescence**: light is absorbed by the molecule and emitted later. Ignoring the picosecond delay in fluorescent small molecules (cf. fluorescent lifetime microscopy), this is exclusively a properties of lanthanide elements (e.g. Europium). Lathanide binding domains have been designed, but do no exist in nature.

# Answer
> 2 votes
The green protein is called GFP of 238 amino acid residues (26.9 kDa) and produced by:

1. the jellyfish Aequorea victoria, It has a major excitation peak at a wavelength of 395 nm.
2. the sea pansy, Renilla reniformis of 498 nm.

However, the protein is widespread and not restricted to these species.

# Answer
> 1 votes
There appears to be a dinoflagellate with Luciferase and the diversity of their genes was studied in the Baltic, found described in an article here: https://www.tandfonline.com/doi/full/10.1080/09670262.2016.1160441

This article by Smithsonian: https://ocean.si.edu/ocean-life/fish/bioluminescence describes a lot of Luciferin metabolizing species in the oceans, saying that different species independently evolved the forms they use. It describes some of the systems like endosymbiotic bacteria, special tissue, and photo-proteins.

---
Tags: phylogenetics, phylogeny, gene, blast, systems-biology
---


# Question
Title: Prevent Snakemake from changing wildcards: sample for every rule
EDIT: I posted this as a bug, perhaps wait for the result of that? github

All,

I am using Snakemake to build a pipeline. To define input and output names, I constructed them based on a filename base (from the samples table), a read indicator (like R1, R2) and the fastq file extension, like so:

```
samples = pd.read_csv('samples.tsv', sep='\t').set_index('samples', drop=False)
reads = ['R1', 'R2']
fastq_file_extension = _001.fastq.gz

```

```
rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

```

However, during pipeline execution already I can see strange things happening, for this rule it reports:

```
[Fri Feb 21 12:39:57 2020]
rule fastqc:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz
output: qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.html, qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.zip
jobid: 10
wildcards: sample=0053_P2017BB3S20R_S2_R, read=1

```

As you can see, Snakemake pulled the R part of the read to the sample and set the read to 1 (should be R1). But this still works because it produces output files as Snakemake expects.

But the next rule, trim\_galore\_pe:

```
rule trim_galore_pe:
    input:
        [os.path.join(config['Directories']['fastq_raw_dir'], '{sample}' + read + fastq_file_extension) for read in reads]
    output:
        fastq_read1 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.1_val_1.fq.gz'),
        fastq_read2 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.2_val_2.fq.gz'),
        report_read1 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.1.fastq.gz_trimming_report.txt'),
        report_read2 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.2.fastq.gz_trimming_report.txt')
    params:
        extra="--retain_unpaired"
    wrapper:
        "0.49.0/bio/trim_galore/pe"

```

This goes wrong, during pipeline execution Snakemake reports:

```
[Fri Feb 21 12:40:37 2020]
Error in rule trim_galore_pe:
jobid: 7
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
conda-env: /home/nlv24077/projects/snaqs/snakemake_envs/1b94430b
cluster_jobid: 11904662.osios
Error executing rule trim_galore_pe on cluster (jobid: 7, external: 11904662.osios, jobscript: /home/nlv24077/temp/RNASeqPairedEnd/.snakemake/tmp.ooypuklc/snakejob.trim_galore_pe.7.sh).
For error details see the cluster log and the log files of the involved rule(s).

```

This is understandeble because the file fastq\_trimmed/0053\_P2017BB3S20R\_S2\_.1\_val\_1.fq.gz does not exist, it should be fastq\_trimmed/0053\_P2017BB3S20R\_S2.1\_val\_1.fq.gz, somehow Snakemake added an extra underscore, which you can see in the execution output, under "wildcards" (0053\_P2017BB3S20R\_S2\_).

```
[Fri Feb 21 12:39:57 2020]
rule trim_galore_pe:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz, fastq_raw/0053_P2017BB3S20R_S2_R2_001.fastq.gz
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
jobid: 7
wildcards: sample=0053_P2017BB3S20R_S2_

```

So how can I prevent Snakemake from changing the samples and read wildcard in seemingly arbitrary ways?

# Answer
> 1 votes
Placing two wildcards directly behind each other will make it impossible for snakemake to determine where one wildcard ends and the other begins. So it is simply taking the longest possible match for the first wildcard (`sample`) and then one character for the second (`read`). This will also cause the `_` to be included in the wildcard value `sample`.

I'm assuming that this same pattern is used for the inputs of another rule to refer to the trim\_galore\_pe output. Hence why the `_` is included in the wildcards there.

You might be able to fix this, simply by adding a `_` in-between `{sample}` and `{readgroup}` in your first rule (and any other rules using this same pattern):

```
rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}_{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

```

# Answer
> 1 votes
I think you want to prevent snakemake from interpreting your wildcard strings as regular expressions. For this, escape each wildcard value using `wildcard_constraints`. E.g., towards the top of your Snakefile add:

```
wildcard_constraints:
    read= '|'.join([re.escape(x) for x in reads]),
    sample= '|'.join([re.escape(x) for x in samples])

```

(Personally, I use the `wildcard_constraints` directive as above quite liberally to constrain every wildcard)

---
Tags: snakemake
---


# Question
Title: Database with antigens and current vaccines
I am seeking virological databases focussing on:

* vaccines
* antibody interactions
* host receptors.

Preferably the data base will host multiple strains/samples of a given virus samples.

Any assistance would be welcome.

# Answer
Given "any virus will do", the viruses that are "head and shoulders" above the rest in terms of volumes of data are HIV and influenza.

The influenza database is here

The influenza viral spike that attaches to the cell receptor is the HA protein – hemagglutinin. The cell receptor is sialic acid – a small sugar that is attached to many different proteins on the cell surface. You can of course simply search PDB for HA or else Genbank.
> 1 votes

---
Tags: phylogenetics, phylogeny, database, covid-19, sars-cov-2
---


# Question
Title: How to download FASTA sequences from NCBI using the terminal?
I have following accession numbers of the 10 chromosomes of *Theobroma cacao* genome.

```
NC_030850.1
NC_030851.1
NC_030852.1
NC_030853.1
NC_030854.1
NC_030855.1
NC_030856.1
NC_030857.1
NC_030858.1
NC_030859.1

```

I need to download these FASTA files using the terminal because I'm working on a server. I'll be grateful if someone can help me. Thank you.

# Answer
Alternatively, you can use the **NCBI Entrez Direct UNIX E-utilities**

Basically, you have to download the install file here: https://www.ncbi.nlm.nih.gov/books/NBK179288/bin/install-edirect.sh

In the terminal, install it using: `source ./install-edirect.sh`

Then, you can download your sequence by doing:

`esearch -db nucleotide -query "NC_030850.1" | efetch -format fasta > NC_030850.1.fasta`

And you should find your fasta sequence downloaded.

As you have several sequences to download, I think it will be quite easy to add this command into a little bash script to process all of them.

More readings about that: http://bioinformatics.cvr.ac.uk/blog/ncbi-entrez-direct-unix-e-utilities/ and https://www.ncbi.nlm.nih.gov/books/NBK179288/
> 11 votes

# Answer
Just `wget` or `curl` each as `https://www.ncbi.nlm.nih.gov/nuccore/NC_030853.1?report=fasta&log$=seqview&format=text`...
> 4 votes

# Answer
## Downloading entire genomes

The best way to download FASTA sequences for an entire genome is to search for the genome, for example `Theobroma cacao genome` in the NCBI Assembly portal and use the big blue Download button.

## Downloading individual chromosomes

For a given assembly, if you want to download the FASTA sequences for a bunch of chromosomes, you can do that by going to the Genomes FTP path highlighted in the screenshot:  Once you are in the Genomes FTP path, you can navigate to the FASTA folder as follows: Theobroma FTP directory \>\> Assembly Structure \>\> Primary Assembly \>\> Assembled Chromosomes \>\> FASTA.

## Downloading a few sequences

For this, you can use Entrez Direct as mentioned by @dc37. Why not always use Entrez Direct? While it is fine for a small number of sequences, it can be slow to download a large number of sequences. Entrez Direct by default will download uncompressed data so you will end up spending more time downloading a larger file instead of downloading a smaller, compressed file from FTP more quickly.

If you were to use Entrez Direct for this purpose, I'd not bother with a bash script and use `epost` to first post the entire list of accessions and then pipe it to `efetch` as shown below:

```
epost -db nuccore -input <file_with_accs> -format acc \
  | efetch -format fasta > sequences.fasta

```
> 2 votes

# Answer
1. Go to the main FASTA sequence page, view source, and get the "ncbi\_uidlist" value. In this case, it's `1061361601`
2. Use the ID in the query string to download directly:

https://www.ncbi.nlm.nih.gov/sviewer/viewer.fcgi?id=**1061361601**&db=nuccore&report=fasta&retmode=text&withmarkup=on&tool=portal&log$=seqview&maxdownloadsize=100000000

You can also get this link directly, by using Chrome's developer tools (F12), viewing the Network tab, then loading the page in (1).
> 2 votes

# Answer
(Sorry, I can't post a comment with my reputation score). @vkkodali, interesting command I did not know about it ;). However, your command is downloading all sequences from the input file into a single fasta file. I believe that if you have a big bunch of sequences, it could be a little bit tricky after to manipulate that kind of file. Just by curiosity, do you have a trick to create one sequence per file using `efetch` ?
> 0 votes

---
Tags: fasta, phylogenetics, linux
---


# Question
Title: Distinguishing Splicings that Produce Functional Transcripts Amongst Alternates Shown
This is a link to the UCSC Genome Browser showing a portion of Ch 17 containing TP53 here.

There are several representations of the gene of varying length and number of exons. I presume these are equivalent to alternate splicings.

Assuming this presumption is correct(?), how can I determine which result in RNAs that are experimentally verified?

# Answer
There are two resources worth checking: GTEx portal and Ensembl. Uniprot and NCBI will call "canonical" the longest sequence comprising valid exons. These are not necessarily the biologically valid ones.

According to Ensembl (TP53, there are a bunch that are experimentally validated and are all class 1. Note the Uniprot column being nearly homologeous: to match the transcripts one has to look at the length, in this case TP53-201 is 393aa long, which is class 1. As mentioned this is not always the case.

So which are more expressed in tissue X? In GTEx, under the collapsed div labelled "Exon Expression" press the "isoform Expression" tab. The most expressed is ENST00000620739.4 (354aa long) join with ENST00000269305.8 (393aa).

In the Exon Expression you will find out what exons are the issue.

Now. to get the sequence of a transcript or covert a variant Ensembl drops the ball. There is no clear web tab that I know of. Nor an API route. ftp://ftp.ensembl.org/pub/release-99/fasta/homo\_sapiens/pep/Homo\_sapiens.GRCh38.pep.all.fa.gz contains all the sequences however. For a project of mine (`ENSTMapper` in this link) to map an amino acid mutation on one transcript to another I ended up write my own short converter, although there are a few similar tools out there.
> 0 votes

---
Tags: genome-browser
---


# Question
Title: tool to visualize a collection of phylogenetic trees
I have a collection of phylogenetic trees produced from clustering DNA sequences, mostly a few nodes in each cluster/tree, with several clusters of size 1. And I would like to visualize them all together to get an idea of the relative sizes of the trees in the collection of clusters (cluster set).

E.g.

```
(ABC:0.0,('CL-118070':0.02010925,'CL-117615':0.00526463):0.00443474);
(DEF:0.0,('CL-117360':0.00967856,'CL-117181':0.00967856):0.00990975);
(GHI:0.0,('CL-123443':0.00123443,'CL-123443':0.00132446):0.00123435);
(JKL:0.0,('CL-134123':0.00123441,'CL-545443':0.00343446):0.00134335);
((((CL-114948:0.00908659,(CL-114989:0.00324875,CL-114959:0.00536289):0.00037038):0.00457230,(CL-114851:0.01291841,CL-114831:0.00590321):0.00077034):0.03411459,XYZ:0.00251364):0.00207594,(CL-114950:0.03074481,CL-114699:0.02256162):0.00733986,CL-116426:0.01643848);

```

I am thinking there must be a tool where I give it an input file with a tree per line and shows them as all in one go.

Any suggestions:

EDIT: looking at the options below:

https://en.wikipedia.org/wiki/List\_of\_phylogenetic\_tree\_visualization\_software

FigTree: in this tool, one can load several trees, one per line, and there is a 'Prev/Next' menu to iterate over them, but no apparent option to display all of them at once.

Phylo.io: in this tool, one can load several trees, and display one at a time or two of them, with colouring options that compare the common nodes/branches. There isn't an apparent option to visualize many trees in one go.

PhyBin: this tool (https://cs.indiana.edu/~rrnewton/projects/phybin/) seems to cluster the sets of trees, but not visualize all of them together.

GrapeTree: this tool () can visualize only one tree at a time. It displays the last tree in the list if more than one are copy+pasted in the input text box.

There is another group of tools under Graph drawing: https://en.wikipedia.org/wiki/Graph\_drawing

Cytoscape: this tool is designed for network visualizations, which could be an option. Unsure if it visualizes newick trees, but there is are tools to convert trees to networks that could potentially work in combination with Cytoscape, e.g. https://github.com/MonsantoCo/Newick2Network

# Answer
I recommend using `densitree`: https://www.cs.auckland.ac.nz/~remco/DensiTree/

And `rwty`: https://github.com/danlwarren/RWTY

Densitree plots the superimposed set of trees in a heatmap-like format:

And `rwty`, a very powerful `R` package that allows visualisation of treespace in an NMDS format (among other methods):
> 2 votes

# Answer
One possible solution is to use ETE Toolkit Treeview, which given a set of Newick files, it will show all of them, albeit connected as shown in the screenshot:
> 1 votes

# Answer
For comparing a large number of trees side-by-side - Dendroscope.

For combining the trees into a network - Splitstree.
> 0 votes

# Answer
Dendroscope is worth a look. I would also look into the R phangorn package.
> -1 votes

---
Tags: phylogenetics, visualization, phylogeny
---


# Question
Title: Program Quality control for data sequencing (nanopore)
Hope all of you are fine. I am wondering if exists a program easy to use that can bring me Basic Statistics of Quality Control from raw data made by MinION device (nanopore).

I have read about FASTAQC and Porechop. My sequencing raw data was a result of a metagenomics analysis (16s and 18s). My main objective is only to filter the sequences and obtain a basic statistic about the lenght and the number of sequences ordered by the lenght. Thank you ladies and gentlemen

PD Sorry for my poor english. Day to day i'm trying to improve it.

# Answer
A simple start would be to use MinIONQC, which is also available via bioconda under the name `r-minionqc`.
> 1 votes

# Answer
I assume raw data is in Fastq format, but ideally you should also use the sequencing\_summary.txt file generated by the basecaller. Other alternatives for QC are PycoQC and NanoPlot. I developed and maintain NanoPlot, so let me know if there is something I can help you with.
> 1 votes

---
Tags: nanopore
---


# Question
Title: Trying to divide number in fourth TSV column by a certain number according to its corresponding identifier in the first column
I was wondering if anyone could help with my little issue in manipulating tsv files in order to calculate estimated coverage across many tsv files.

I have tsv files that look something like this:

```
ID        Info    Info2  x
M80428    info    info   1000
T83241    info    info   8400
E81292    info    info   10200

```

I also have another tsv file with all possible IDs and a corresponding number like so:

```
ID        y
M80428    27311
E11123    46531
P38163    14213
A34197    12892
T83241    37416
Q21345    27321
E81292    13101

```

x in the first tsv are estimated read counts and the y in the second is the genome length of each ID.

I would like to divide the "x" by the "y" in the second tsv that has the corresponding ID in order to get the estimated coverage of each hit. e.g. the first row would be 1000/27311, second row 8400/37416. I would then want this to be written to the fifth column of the tsv.

I am working in bash shell mostly here, and think that the way I would probably do this is with awk and by making the second tsv an associative array but I can't work out how exactly to do so. The python package dict seems like another option but my Python skills are practically zero.

Thanks for the help in advance!

# Answer
> 3 votes
Python dictionaries are indeed a good tool for that.

Here is a script that takes 2 arguments (the file containing the counts and the file containing the genome length info), and prints the input data, with an added column containing coverage:

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

counts_fname = sys.argv[1]
genome_lengths_fname = sys.argv[2]

genome_length_dict = {}
with open(genome_lengths_fname) as genome_lengths_f:
    # Checking the header line, for extra safety
    header = genome_lengths_f.readline().strip().split("\t")
    assert header[0] == "ID"
    assert header[1] == "y"
    # Reading the data lines
    for line in genome_lengths_f:
        [genome_id, genome_length] = line.strip().split("\t")
        genome_length_dict[genome_id] = int(genome_length)

with open(counts_fname) as counts_f:
    print("ID\tInfo\tInfo2\tx\tx/y")
    # Checking the header line, for extra safety
    header = counts_f.readline().strip().split("\t")
    assert header[0] == "ID"
    assert header[1] == "Info"
    assert header[2] == "Info2"
    assert header[3] == "x"
    # Reading the data lines
    for line in counts_f:
        [genome_id, info1, info2, counts] = line.strip().split("\t")
        # /!\ It would be different in python 2:
        # different print and division behaviours
        print(
            genome_id, info1, info2, int(counts),
            int(counts) / genome_length_dict[genome_id],
            sep="\t")

sys.exit(0)

```

Testing it:

```
$ ./coverage.py x.tsv y.tsv
ID  Info    Info2   x   x/y
M80428  info    info    1000    0.0366152832192157
T83241  info    info    8400    0.22450288646568314
E81292  info    info    10200   0.7785665216395695

```

# Answer
> 1 votes
I would do it in R:

```
tab1 <- read.table('file1.tsv', header = T, stringsAsFactors = F)
tab2 <- read.table('file2.tsv', header = T, stringsAsFactors = F)
merged_tav <- merge(reads_tab, genome_len_tab)
merged_tav$z <- merged_tav$x / merged_tav$y
merged_tav

```

returns

```
      ID Info Info2     x     y          z
1 E81292 info  info 10200 13101 0.77856652
2 M80428 info  info  1000 27311 0.03661528
3 T83241 info  info  8400 37416 0.22450289

```

note that only IDs that are in both tables will be preserved.

---
Tags: python, bash, awk
---


# Question
Title: How to download multiple proteomes at once?
I'm looking for a way to download multiple proteomes at once from one clade, as fasta. Probably from the NCBI, because it looks more user-friendly than Ensembl, but Ensembl is okay too. In the best case scenario, I would like to constrain the download to a specific clade (eg. Actinobacteria, Chloroflexi etc.), maybe with a specific RefSeq Category (eg. reference/representative genome) or even depending on the genome assembly level. Command line would probably be the best option, but if there's a direct download option it would be nice to know it as well. Does anything like this exist? Help would be much appreciated!

# Answer
You could try download-refseq-genomes. It will fetch all genomes from the NCBI FTP server that are in a specific subtree of the phylogeny.

For example, downloading the amino acid sequences from all proteins in each genome assembly for all species belonging to Chloroflexi (NCBI taxon ID = 200795):

```
download-refseq-genomes.pl -t faa 200795

```

This will download all assemblies with status=Complete. Use option `-a` to download all assemblies with any status.
> 1 votes

# Answer
You can use NCBI resources to download these data directly. For the entire proteome, you can download the sequences of all proteins annotated on the genome assembly in FASTA format. In the NCBI FTP directories, these files have the names with the suffix `_protein.faa.gz`. More details about the files in the FTP directories are described in the README.

## Command-line method using EDirect

This method uses Entrez Direct to first query the Assembly database for all available genomes, parses the output to extract the ftp path and downloads data using wget command.

```
esearch -db assembly -query 'txid200795[organism] AND lastest_refseq[filter]' \
  | esummary \
  | xtract -pattern DocumentSummary -element FtpPath_RefSeq \
  | while read -r url ; do
    path=$(echo $url | perl -pe 's/(GC[FA]_\d+.*)/\1\/\1_protein.faa.gz/g') ;
    wget -q --show-progress "$path" -P protein_data ;
  done

```

Here, I am restricting the results to latest RefSeq data by using the `latest_refseq[filter]` but you can use any other Entrez filters to further narrow down your results. Along the same lines, you can modify the code minimally as shown below to download multiple file types:

```
esearch -db assembly -query 'txid200795[organism] AND latest_refseq[filter]' \
  | esummary \
  | xtract -pattern DocumentSummary -element FtpPath_RefSeq \
  | while read -r url ; do
    base_path=$(echo $url | perl -pe 's/(GCF_\d+.*)/\1\/\1_/g') ;
    for suf in genomic.fna.gz protein.faa.gz ; do
      wget -q --show-progress "$base_path$suf" -P ncbi_data ;
    done ;
  done

```

## GUI method

You can use the same query as shown above in the NCBI Assembly portal and then use the 'Download Assemblies' button to pick and choose the type of files you would like to download. By default, if you do not select any assemblies in the results, this interface downloads files for all results.
> 0 votes

---
Tags: phylogenetics, ensembl, data-download, proteomics
---


# Question
Title: Importance of Proper Pairs vs Aligned Reads for RNASeq data
I have stranded, paired end RNASeq reads that I have aligned using STAR. I plan do conduct a differential expression analysis with DESeq2. After running quality control checks, a good portion of my reads are aligned but I'm a little concerned about the proportion of properly paired reads

Here is the inner distance plot for the samples

For human RNA sequencing, does anyone have a rule of thumb that they'd consider is an appropriate proportion of properly paired reads? Are there any good papers that talk about this issue? As MultiQC marks the column as yellow, I'm concerned the proper pair percentage may not be satisfactory.

How important are properly paired reads for differential expression analyses?

# Answer
> 1 votes
Differential gene analysis sees only the read counts. So the proportion of not properly paired reads will have an effect on the counts, if the programme / software used, counts fragments / read pairs instead of single reads. For example HTSeq will count

> From https://htseq.readthedocs.io/en/release\_0.11.1/count.html
>
> For paired-end data, does htseq-count count reads or read pairs? Read pairs. The script is designed to count “units of evidence” for gene expression. If both mates map to the same gene, this still only shows that one cDNA fragment originated from that gene. Hence, it should be counted only once.

In Rsubread, featureCounts you set `requireBothEndsMapped=TRUE` to count fragments.

Short anaswer is if the sequencing is deep enough, you should still have enough reads. In your situation, I estimated that you lose about 40% of your reads if you use only properly paired reads, which means you lose quite a bit of counts.

The larger question is what reads are you missing?

I see that you have about 20-30% of reads that are non-uniquely mapped, I suppose part of this will contribute towards the not properly paired reads. I would check first whether most of the not properly paired comes from multiple aligned.

Then if not, you can check where these not properly paired reads are falling, if they are some repetitive regions or say rRNA etc, then most likely thats a problem with the library prep..

# Answer
> 2 votes
If you check the RSeQC docs you can get an explanation of how the inner distance stats work. By the look of your inner-distance MultiQC plot I guess that your sequencing is 2x100bp reads. This plot shows that your insert sizes on the pairs is pretty short, meaning that a lot of the paired reads are overlapping significantly. This will be affecting the alignments a little, resulting in the low proportion of properly paired reads.

The good news is that for differential expression, this probably doesn't really matter too much. What is important is getting fragment counts, so as long as at least one of the pair mates is mapping properly, then you will will be getting sensible counts - even if read 2 is not mapping properly. You should be able to have confidence in your results from multiple replicates if in doubt - you have more than just `-rep1` right? ;)

Note that the yellow background in the MultiQC column is not passing judgement with categories like pass/warn/fail in FastQC - it just happens that the colour scheme ended up there for your samples. It's a continuous gradient from red to yellow to green for those values (0-100%).

---
Tags: rna-seq, ngs, transcriptome, quality-control, star
---


# Question
Title: vcf file to gene annotations and or genes
Is there some tool that will convert a vcf/bcf file to some gene annotation, based on the fact that I have the gene annotation of the already exsisting reference genome. Or some tool that will give me the genes based on the most probable chance.

# Answer
Creating a FASTA file by applying variants found in a VCF file is called "creating a consensus FASTA file". You can use `bcftools consensus`. Unfortunately, the tools only works genome-wide (at the FASTA level), so you might need to use a different tool to extract a gene sequence from the FASTA file using the gene's co-ordinates.

If the VCF file is annotated with HGVS notations for the variants, there might be a way to use those notations to apply the edit to the reference gene/transcript sequence directly. I know the hgvs python module works well with the HGVS notation, but I don't see a way to "apply" the variant to the reference sequence.
> 2 votes

---
Tags: vcf, gene
---


title: Joining columns to a sorted file

text:
# Question
Title: Joining columns to a sorted file
I have had a segmentation file (copy number) like below

```
  Sample      Chromosome  Start      End          Segment_Mean
LP6008334-DNA_C02   1     10583      62271      -0.992452416
LP6008334-DNA_C02   1   17201069    17230953    -0.992452416
LP6008334-DNA_C02   1   17235034    17235196     0.592510084
LP6008334-DNA_C02   1   142781738   142967188   -0.992452416

```

I was using bedtools to count the number of markers in above genomic range (Start, End) when I got as below

```
1   10583   62271   6
1   10583   62271   6
1   10583   750153  75
1   14907   1062025 106
1   121009  245455  13
1   404187  5719467 532
1   732809  5726842 500
1   821926  1120377 31
1   952428  1179983 23
1   1181751 1469134 29
1   3458681 143542402   14010
1   5726844 5736229 2
1   5735836 16832446    1111
1   5736682 8081290 236
1   8353020 8476455 13

```

Basically the same files but the second one has extra last column for the number of marker but sorted one I then I can not matched this with the above sample names and Segment\_Mean

How I can add the corresponding sample names and Segment\_Mean to my second file?

I have tried this where I am getting this mess

```
DF <- left_join(df_res1,df, by = "Start")

> head(DF)
           Sample.x Chromosome.x     Start     End.x Num_Probes          Sample.y Chromosome.y     End.y
1 LP6008334-DNA_C02            1     10583     62271          9 LP6008334-DNA_C02            1     62271
2 LP6008334-DNA_C02            1     10583     62271          9 LP6008337-DNA_A07            1     62271
3 LP6008334-DNA_C02            1     10583     62271          9 LP6008460-DNA_F02            1    750153
4 LP6008334-DNA_C02            1  17201069  17230953         20 LP6008334-DNA_C02            1  17230953
5 LP6008334-DNA_C02            1  17235034  17235196          1 LP6008334-DNA_C02            1  17235196
6 LP6008334-DNA_C02            1 142781738 142967188         40 LP6008334-DNA_C02            1 142967188
>

```

Also with this I have some overlaps like

```
300 segment overlaps detected in file '/temp/hgig/fi1d18/man/examplefiles/res1.txt'.

```

First overlap detected between segments at lines 1 and 9.

# Answer
> 2 votes
Rearrange your text file into a sorted BED file:

```
$ tail -n+2 segmentation.txt | awk -vFS="\t" -vOFS="\t" '{ print $2, $3, $4, $1 }' | sort-bed - > segmentation.bed

```

Use BEDOPS `bedmap` with `--echo-map-id-uniq` and `--count` operators on merged regions:

```
$ bedmap --delim '\t' --echo --count --echo-map-id-uniq --fraction-map 1 <(bedops --merge segmentation.bed) segmentation.bed > answer.bed

```

The file `answer.bed` will contain merged intervals from `segmentation.bed` in the first three columns, the number of intervals from `segmentation.bed` which are contained entirely within the merged interval in the fourth column, and a unique listing of IDs of those contained intervals in the fifth column.

---
Tags: r, bedtools, wgs, cnv
---


# Question
Title: Why does samtools mpileup sometimes include ref bases (other than ',' or '.')?
This is my first post here. I can think of no way of giving an easily reproducible example, as per the stack ethos. SO apologies in advance and any feedback on question format appreciated.

I have called samtools mpileup (Version: 1.9 (using htslib 1.9)) on some oxford nanopore data (mapped to GRCh38).

I am seeing some instances like below (grouped with pythons Counter), here reference is T at this position:

```
 Counter({',': 34, '.': 33, 'c': 9, 'g': 6, 'a': 4, 't': 2, '-:ACTC': 1, '-:ACTCGG': 1, '-:5': 1, '$': 1, '*': 1})

```

Why is 't': 2' in there?

flags used: -A -B -Q 1 -R -x

# Answer
> Why is 't': 2' in there?

How else do you think the pileup should show you that two reads (in the reverse direction) have a T in that position?
> -1 votes

---
Tags: python, samtools, long-reads, mpileup
---


# Question
Title: Are there any tools for discovering tandem repeats around a variant in a VCF file?
I know for example `PolyX` from Pierre's `Jvarkit` to find the number of repeated REF bases (homopolymers) around a position in a VCF file as described here exists.

I would like to know however if there tools that describes tandem repeats around a reported variant (TR) up to 6 repeated bases. I can imagine 3 different scenarios:

**Scenario 1:**

```
chrK    posY    C     G

TR of (ATC)
____________v
ATCATCATCATCGATGATAC

```

**Scenario 2:**

```
chrK    posY    A     T

             TR of (AACT)
            v____________
CAGTCGATTGGATAACTAACTAACT

```

**Scenario 3:**

```
chrK    posY    C     A

TR of (AG)    TR of (AG)
____________v__________
AGAGAGAGAGAGCAGAGAGAGAG

```

If you know any tools in mind that can be used for VCFs, I'd really appreciate citing them.

Thanks in advance.

# Answer
> 1 votes
After researching a bit, I found a couple of tools that take as input BAM files. I did not test them yet:

If I find more tools, I will update this post. After testing them, I will come back to share my feedback.

# Answer
> 1 votes
The old program for doing this last updated in 2016 is Tandem Repeat Finder (TRF) here. This will not find a microsatellite sequence at a given variant, but instead find all tandem repeats in the genome and their position. It has been around for a long time.

You would then simply screen this output against the position of your variant. TRF may have been modified however to perform a search within given nucleotide boundaries.

---
Tags: vcf, repeat-elements
---


# Question
Title: K means clustering, would PCA be a better option?
I have the data below. I need to use a clustering method to classify them and into categories of "Heterozygotote, Allele 1, Allele 2 and No Call. The values in RFU1 and RFU2 are used to determine the call variable.

Any values in RFU1 and RFU2 that are negative should automatically be classified as No Call. I want to develop my own validation check to check the PCR program calls.

**Data**

```
structure(list(Well = structure(1:96, .Label = c("A01", "A02",
"A03", "A04", "A05", "A06", "A07", "A08", "A09", "A10", "A11",
"A12", "B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08",
"B09", "B10", "B11", "B12", "C01", "C02", "C03", "C04", "C05",
"C06", "C07", "C08", "C09", "C10", "C11", "C12", "D01", "D02",
"D03", "D04", "D05", "D06", "D07", "D08", "D09", "D10", "D11",
"D12", "E01", "E02", "E03", "E04", "E05", "E06", "E07", "E08",
"E09", "E10", "E11", "E12", "F01", "F02", "F03", "F04", "F05",
"F06", "F07", "F08", "F09", "F10", "F11", "F12", "G01", "G02",
"G03", "G04", "G05", "G06", "G07", "G08", "G09", "G10", "G11",
"G12", "H01", "H02", "H03", "H04", "H05", "H06", "H07", "H08",
"H09", "H10", "H11", "H12"), class = "factor"), Sample = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
    Call = structure(c(4L, 4L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L,
    3L, 3L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 4L,
    4L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 1L, 1L,
    3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L,
    1L, 2L, 2L, 1L, 1L, 3L, 3L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
    2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L), .Label = c("Allele 1",
    "Allele 2", "Heterozygote", "No Call"), class = "factor"),
    Type = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "Auto", class = "factor"),
    RFU1 = c(-0.295502405, 0.964070798, 3381.332182, 3532.769062,
    3431.836843, 3242.966511, 2104.791167, 2220.008503, 3548.252161,
    3506.51418, 2290.273178, 2281.587684, -5.64819475, -11.73109864,
    3784.914039, 3619.00781, 3618.211608, 3248.106466, 3394.650325,
    3339.870196, 2449.202902, 2426.835174, 3432.153478, 2761.392304,
    -9.267907504, -7.365704356, 3743.092314, 3787.241702, 2172.027787,
    2096.845649, 2135.649551, 2149.145547, 2293.757257, 2348.099108,
    2321.019045, 2022.168867, -17.93532331, -12.59832941, 3805.416768,
    3498.998138, 2304.597239, 2509.63987, 2181.11547, 2261.011876,
    3432.453036, 3662.758933, 2371.11049, 3068.827061, 2612.107589,
    2687.824075, 3179.315918, 3688.525218, 3465.327523, 3405.154043,
    2535.514915, 2452.200079, 374.435299, 423.6015308, 3742.515563,
    3578.777925, 2634.955017, 2527.514043, 3817.579252, 3550.999412,
    -10.72035816, 3294.486334, 3352.40368, 3463.150507, 3472.576514,
    3741.898759, 3571.369947, 3720.645869, 3739.569593, 3855.583168,
    418.6837047, 49.47548241, 2171.034284, 2155.314713, 3432.363384,
    3582.508917, 3425.415274, 3487.203299, 3505.23909, 3413.342096,
    113.5100691, 128.6414053, 2454.588175, 2323.061591, 3188.705702,
    3376.950546, 3291.072437, 3181.001961, 3195.013863, 3776.919197,
    2284.22659, 2277.338631), RFU2 = c(-8.346468029, 235.4058561,
    637.9218251, 650.3759507, 617.4161748, 604.0792911, 4270.310727,
    4199.615749, 689.863543, 712.6144338, 4274.287194, 4541.168491,
    -1.626221758, -2.437395631, 802.0941252, 730.5998997, 686.9037384,
    625.8245403, 644.3586836, 642.8833044, 4937.691887, 5159.479928,
    725.4449756, 573.3910899, -4.006398006, 213.2859144, 739.7910786,
    731.0150586, 4408.81923, 4767.533969, 4302.641493, 4325.913445,
    4597.47663, 4666.904418, 4800.357526, 4142.535329, -17.23239968,
    178.5311942, 778.305843, 743.1438168, 4214.507094, 4553.703511,
    4629.339014, 4459.697405, 661.7299014, 727.1054982, 4553.170272,
    5482.231486, 4520.517999, 4737.802036, 643.3599887, 726.4314715,
    696.5968338, 697.6099599, 411.8118071, 409.4943424, 5687.32635,
    5757.51512, 766.4240193, 779.2403225, 4745.055632, 4582.267792,
    749.5679421, 675.8747055, -7.254521898, 628.3467565, 631.116767,
    672.7064514, 687.2642132, 718.1192917, 731.785499, 668.3686048,
    784.8055727, 791.3155894, 4471.047168, 4501.597841, 4504.670332,
    4442.621066, 682.0632225, 706.6204595, 680.5242182, 683.9558692,
    684.2909706, 618.6535251, 5727.684954, 6098.485474, 5099.952926,
    4779.742057, 571.4303822, 614.9258218, 602.9830491, 651.2847695,
    591.8833499, 742.2387568, 4443.376841, 4716.792177)), class = "data.frame", row.names = c(NA,
-96L))

```

**What I have tried so far**

```
library(cluster)
library(factoextra)
library(formattable)

df <- df[,c(1,5,6)]
df$RFU1[df$RFU1 < 0] <- 0
df$RFU2[df$RFU2 < 0] <- 0
df$RFU1 <- formattable(df$RFU1, digits = 2, format = "f")
df$RFU2 <- formattable(df$RFU2, digits = 2, format = "f")
df$Well <- as.numeric(df$Well)

clusters <- kmeans(df, centers = 4)
Kmeans_plot <- fviz_cluster(clusters, data = df)

```

**This is the plot generated**

Points in the top right 57,58,75,76,85,86 should be in another cluster (classed as allele 2 but placed in cluster 1 (rest of points in this cluster are heterozygote)

Also points 24, 55, 56 should be clustered with cluster 2 as all these points are classed as allele 1

**I need an algorithm that can detect optimal cluster center, some PCR runs may only have 3 different call variables instead of 4**

**Programs output - Software company stated their algorithm isn't sophisticated**

# Answer
I would stay away from using k-means, and instead use a method that doesn't *a priori* define a number of clusters to detect. It also looks as though your clusters aren't exactly spherical, which is an assumption of k-means. I personally am a fan of `dbscan`, which is available in the R package of the same name. The other poster recommended t-SNE (available in the `Rtsne` package) but wasn't exactly clear on how / why to use it. t-SNE is a nonlinear dimension reduction technique used specifically for visualization, so you could cluster your points in *n*-dimensional PCA space and then plot the results in two-dimensional t-SNE space. I think the other poster was recommending plotting your data in t-SNE space before clustering to get an idea of how many true clusters are available, which is another option you could pursue.
> 1 votes

# Answer
What happens when

```
clusters <- kmeans(df, centers = 3)

```

Anyway K means is pretty good and hasn't done a bad job of your data, although I don't know what you mean by the "cluster centre". There are methods to establish 'centroids'.

With K-means you need to vary the group size or have a very good a priori why it is 4, rather than any other number.

The analysis to perform is tSNE ... either by itself or else following PCA. This will provide a clear a priori of the number of groups in your analysis. I don't know R, but it is present in scikit-learn in Python.

---

Your observation that some of points are in the wrong cluster ... increase the number of groups in k-means, or switch to PCA, tSNE or PCA-tSNE. If the points remain in the "wrong groups" in all analyses at different group sizes for k-means, then basically thats the result. K-means is good, but not perfect, and misclassification will happen.

Another route is the ML repertoire of lasso regression, ridge regression, decision trees, random forests etc ... they all do the same basic job but I'd look at K-means, PCA, tSNE first and if ok ... last. Regression style calculations will give different answers.

---

t-SNE is a dimensionality reduction approach to produce a bivariate plot and the fondness of this approach is purely observational rather than empirical. It works and also works very well following PCA. Its the method used in Seurat. K-means is used in ML, and certainly has a fan base, but it is not to say it will give you a correct answer.
> 0 votes

# Answer
I don't think clustering algorithms are working here, it looks like you can do this much simpler by looking at whatever determines the X axis in your plot. The samples split cleanly when looking in just one dimension.
> 0 votes

---
Tags: r, clustering, k-mer, pca, pcr
---


# Question
Title: Identifying mutually **exclusive** gene sets
**I am interested in identifying gene pairs (or better: sets of genes) whose expression is mutually exclusive.** Ideally, both genes (or gene sets) would be widely expressed but I am also interested in cases that are very cell type specific.

Familiar examples of the former would include some of the cyclin genes, in humans e.g. cyclin E and cyfclin B, which are expressed during S-phase and M-phase, respectively. Another example would be some of the clock-genes, where Per1 and Per2 are often highly expressed in the suprachiasmatic nucleus (SCN) before midday, and whereas Bmal1 is expressed after midday. However, such genes need not be cyclically expressed, they can simply be very specific markers for certain cell types. For example, GAD65 and GAD67 are very specific for interneurons, whereas CamKIIa is very specific for excitatory glutamatergic cells.

I have been able to find a single paper so far that systematically tries to find such genes from gene expression data, namely: Tzanis & Vlahavas (2010) Mining for Mutually Exclusive Gene Expressions. According to Google Scholar, this paper is cited by a single other paper that does not seem to be relevant. In Tzanis & Vlahavas (2010), they use a SAGE library (from 2002) of short expression sequence tags to find these strong negative associations. By today's standards, the data set is pretty small (basically, it contains 90 samples). The paper itself is hence anything but a canonical source for such negative interactions.

### Questions:

1. Is there a curated database of mutually exclusive gene pairs or gene sets?
2. If not, is there a curated database that would allow me to find such genes?

Otherwise, my current plan would be to use the human or mouse cell atlas single cell data to define a signed, weighted gene expression network. I would then search for pairs of clusters/modules, where genes have predominantly positive correlations within a cluster but predominantly negative correlations between clusters. So essentially, signed, weighted gene correlation network analysis (WGCNA) with extra steps.

3. Has somebody published such an analysis already (the WGCNA literature is pretty substantial)?
4. Is there other (recent) related work that I am unaware off?

I am fairly new to bioinformatics, so the lack of papers and resources that I am able to find is likely mostly a reflection of my lack of ability to find the correct search term.

### Addendum

Similar biostars thread from 5 months ago, also without an answer.

# Answer
> 1 votes
"Mutually exclusive" is not a precise statistical term because it could be seen as independence.

Anyway if truly mutually exclusive AND their behaviour follows periodicity, i.e. sine waves, which it might do in a cell cycle, i.e. the cells divide at periodic intervals, you can do this via asynchronicity.

If the gene expression can be described under periodicity analysis you simply sum the amplitude between two genes. When the amplitude is zero ... they are expression is antogonistic - when one is expressed the other is silent.

Amplitude simply means the height of the period.

Periodicity analysis is quite simple if you have uniform time measurements.

If you have 'missing data', i.e. you time intervals are not uniform, for example you skipped a reading then this is called missing data in periodicity analysis. Missing data periodicity analysis is complicated, but its doable.

So in this context its quite simple, you would identify the periodicity of all genes and construct a 2x2 matrix subtracting their respective amplitudes. It would be kinda juicy to do.

---

The above is cool if the genes are periodic, however covariance is a better solution in general terms.

What you want is **negative co-variance**. If two genes had a covariance of -1 ... that will mutually exclusive. As above 2x2 matrix and the co-variance across the grid. You do need to be careful that the program you are using doesn't remove negative values of the covariance.

```
Covariance = 1, meaning two genes are co-expressed
Covariance = 0, random distribution
Covariance = -1 mutually exclusive expression

```

---
Tags: scrnaseq, networks, wgcna, correlation
---


# Question
Title: Drawing 2 heatmaps with the same gene order
I have 2 dataframe with same genes both different data sets

like

```
> head(heat1[1:3,1:4])
             A10         A11        A12        A2
ATP5F1 2.1163092  0.32044094 -0.5661269 0.7978236
DDX5   0.6277182 -1.67055853 -0.0123082 0.8895472
EEF1G  1.6459212  0.02923401 -1.2737676 0.5842162
>

> head(heat2[1:3,1:4])
              A10        A12          A3          A4
A2M    -0.6031969  0.6378864 0.000573381 -0.07769313
ABCB11  0.5600164  0.3087120 0.171636836  1.25681512
ABCC2   0.2093133 -0.5608016 0.782278709  0.16618683
>

```

I want to plot them side by side with the same order of genes

I tried

```
> ht1 = Heatmap(heat1)
> ht2 = Heatmap(heat2)
> ht_list = ht1 + ht2
> draw(ht_list, ht_gap = unit(1, "cm"))
Warning message:
Row names of heatmap 2 is not consistent as the main
heatmap (1)

```

But I getting not the same order of genes

Can you help me with this?

# Answer
> 3 votes
Please read through the ComplexHeatmap manual, especially the section on row and column orders: https://jokergoo.github.io/ComplexHeatmap-reference/book/a-single-heatmap.html#get-orders-and-dendrograms-from-heatmap

You can use `row_order` and `column_order` following a `draw()` (which you can nest between a `pdf(null)` and `dev.off()` to disable drawing it on any device). That way, you'll be sure to get the right order from the first heatmap to use with the second.

By nesting, I mean like so:

```
h1 <- Heatmap(...)
pdf(NULL)
h1d <- draw(h1)
dev.off()
ro1 <- row_order(h1d)
co1 <- column_order(h1d)

```

Use whichever (row/column) order is appropriate to gene order.

---
Tags: r, complexheatmap
---


# Question
Title: How are these definitions related to differential expression?
I have two groups of patients; for each patient I have an output file (RNA-seq) contains this information

```
                  readcounts_union  readcounts_intersectionNotEmpt  genelength  R/FPKM_union    R/FPKM_intersectionNotEmpt
ENSG00000000003.10             8                      8               3688     1.06120392730222        1.05722142116966

ENSG00000000419.8             31                     31               4032     3.76132572546529               3.74721015127372

```

People say

```
readcounts_union

```

***reads that overlap a single feature are included in the count for that feature; reads that overlap more than one feature are not included in the counts***

```
readcounts_intersectionNotEmpty

```

***regions that are shared between features are excluded; reads that overlap a remaining region are included in the count for the feature that includes that region; reads that overlap regions from more than one feature are not included in the counts***

```
R/FPKM_union

```

***readcounts\_union normalised for gene length and sequencing depth***

```
R/FPKM_intersectionNotEmpty

```

***readcounts\_intersectionNotEmpty normalised for gene length and sequencing depth***

I want to do differential expression analysis between groups but I don't know which column I should use for DESeq2

Anybody familiar with such a format?

# Answer
Either of the first two columns are fine to use, though I personally prefer the "intersectionNotEmpty" one.

As an aside, that looks like a custom format file produced from post-processing htseq-count output.
> 2 votes

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


title: gnu parallel for macs peak calling

text:
# Question
Title: gnu parallel for macs peak calling
These are the list of files

```
Control_Input_sorted.bam_rem.bam
Control_H2BUb_sorted.bam_rem.bam
Control_IgG_sorted.bam_rem.bam
PTPN6_g2_6_H2Bub_sorted.bam_rem.bam
PTPN6_g2_6_Input_sorted.bam_rem.bam

```

So to get unique files i do this

```
 ls -1 *bam | sort | sed -r 's/_sorted.bam_rem.bam//g' | sort | uniq
Control_Input
Control_H2BUb
Control_IgG
PTPN6_g2_6_H2Bub
PTPN6_g2_6_Input

```

Now i have to run each sample against input and IgG for peak-calling. Like `Control_H2BUb_sorted.bam_rem.bam` against

```
Contol_Input_sorted.bam_rem.bam & Control_IgG_sorted.bam_rem.bam and

PTPN6_g2_6_H2Bub_sorted.bam_rem.bam against PTPN6_g2_6_Input_sorted.bam_rem.bam

```

So how do I parse and execute using gnu parallel.

I came across this tutorial where it does something like this

```
cat sample_names.txt | parallel --max-procs=12 'macs2 callpeak -t {}-A-NC.sorted.bam \
 -c {}-G-NC.sorted.bam -g hs -n {}-A-NC-sharp-model -q 0.01 --outdir {}-A-NC-sharp-model-peaks 2> {}-A-NC-sharp-model.stderr'

```

Now im not sure how to pass the file names as argument for IP and Input or IgG.

If i use logic a little bit i have to fix Control and PTPN6\_g2\_6 as constant naming!.

Any suggestion or help would be really appreciated.

# Answer
not gnuparallel but **nextflow**. I don't know macs and what is 'first.bam' and 'second.bam' is not clear not me but I hope you'll get the idea.

The nextflow:

```
Channel.fromPath(params.bams).splitCsv(header: false,sep:'\t',strip:true).map{T->file(T[0])}.into{bams1;bams2}

boolean isControl(file) {
    return file.name.endsWith("_Input_sorted.bam_rem.bam") || file.name.endsWith("_IgG_sorted.bam_rem.bam ");
    }

bams1.filter{T->isControl(T)}.set{controls}
bams2.filter{T->!isControl(T)}.set{experiment}

process mac {
    maxForks 14
    tag "${bam1} vs ${bam2}"
    input:
        set bam1 , bam2  from experiment.combine(controls)
    output:
        file("output.peaks") into mac_peak
    script:
    """
    # check syntax please
    macs2 callpeak -t {bam1.toRealPath()}  -c ${bam2.toRealPath()} -g hs -n out-model -q 0.01 --outdir output.peaks
    """
    }

```

invoke:

```
find /dir1/dir2 -type f -name "*.bam" > input.list

nextflow run --bams input.list workflow.nf

```
> 1 votes

---
Tags: shell, macs2
---


title: Filtering raw sequencing reads

text:
# Question
Title: Filtering raw sequencing reads
I have to filter the raw sequencing reads based on the following criteria:

1. Remove reads containing adapters
2. Remove reads containing N \> 10% (N represents base that could not be determined)
3. Remove reads where the Qscore (Quality value) of over 50% bases of the read is \<= 5

The reads are paired-end fastq files. I have been able to do the first part by using trimmomatic:

```
trimmomatic \
  PE in_1.fq.gz in_2.fq.gz \
  pair_out_1 unpair_out_1 pair_out_2 unpair_out_2 \
  ILLUMINACLIP:adapter.fa:2:30:5

```

The other 2 steps are a bit more challenging since I have not found a tool that perform exactly what is required. I could write some custom code of course but I would like to use the existing tools if possible. No need to reinvent the wheel.

Do you have any suggestions on how to proceed?

# Answer
I think you can come close to this with `cutadapt`:

```
cutadapt \
  -a ${FWD_ADAPTER} -A ${REV_ADAPTER} \         ## adapter
  --max-n 0.1 \                                 ## do not allow > 10% Ns
  -q 5 \                                        ## remove bases with basequal < 5
  --minimum-length (readlength/2) \             ## minimum length after trimming
  -o read1_trimmed.fq -p read2_trimmed.fq \     ## outputs
  read1.fq read2.fq                             ## inputs

```

This will remove reads with more than 10% of bases being `N` and remove bases with base quality below 5. If you then set `--minimum-length` to half the read length, this should also remove reads where more than half of the bases had a base quality below 5 and therefore got trimmed away. Poor base quality in Illumina reads typically happens towards the 3' end so I think this little workaround is a sufficient proxy to achieve your goal. A base quality of 5 is very low in any case, so if sequencing went normal the typical quality should be notably higher. I would be surprised if this third filter is really necessary and will eliminate a notable number of reads.
> 2 votes

---
Tags: fastq, quality-control, filtering
---


# Question
Title: Using cellranger mkgtf w/ NCBI GFF3 files
I am using several cellranger tools/pipelines to analyze some data, specifically, mkgtf, mkref, count, and aggr. I am using mkgtf and mkref to include a PI-synthesized protein to my reference. I have finished the analysis but am not seeing the protein in my output, so I am going back to the start to verify the process.

Unlike the documentation (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references), my GTF file was not downloaded from ENSEMBL or UCSC, but NCBI (https://www.ncbi.nlm.nih.gov/nuccore/AY678269.1?report=genbank). This GFF file does not have the same key nomenclature as seen in ENSEMBL or USCS.

I did use the mkgtf tool, but the created file is identical to the input.

```
cellranger mkgtf custom_ref/AY678269_genes.gtf reference_data/AY678269_mkgtf_genes.gtf

```

I went ahead and used cat to add the fasta and gtf file from the synthesized protein to the respective reference fasta and gene files, then ran mkref.

The s.protein is in the fasta file:

```
>lcl|AY678269.1_cds_AAV52169.1_1 [protein=tandem-dimer red fluorescent protein] [protein_id=AAV52169.1] [location=1..1431] [gbkey=CDS]

```

as well as the gene list

```
##sequence-region AY678269.1 1 1431
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=32630
AY678269.1      Genbank region  1       1431    .       +       .
AY678269.1      Genbank CDS     1       1431    .       +       0

```

However this protein is not found in the output...surely I must have made an error. Does anyone have any ideas about what I could be doing incorrectly?

Do I need to change the entire nomenclature of the GFF keys to match what I see from ENSEMBL (e.g., gene\_id ""; gene\_version ""; transcript\_id, etc...)?

Any assistance is greatly appreciated.

# Answer
> 1 votes
Regarding your comment:

> This GFF file does not have the same key nomenclature as seen in ENSEMBL or USCS

The link you have provided corresponds to an entry in `Genbank` format and not to `GFF` format. I believe you need a file in `GTF` format in particular.

Nevertheless, given the fact that you have manually added the corresponding entries to the reference `fasta` and `gtf` files prior to `mkref`, you should have seen your gene of interest in the output. I would double check one of the original fasta headers and original gtf entries in the reference file and match these for the new entry.

---
Tags: cellranger
---


title: How to run trimmomatic in HPC

text:
# Question
Title: How to run trimmomatic in HPC
Sorry if this has been asked before, I've done a quick search and I don't think there's an easy explanation for me to understand.

I'm really new to bioinformatics/RNASEQ analysis, having only taken up an attempt to process my RNASEQ data last week. I have 0 background and some rudimentary knowledge of command line. So far I've managed to get access to my universities HPC and perform a fastqc on my fastq data.

What I'm stuck at at the moment is attempting to trim my RNASEQ data and remove adapters. My data had been run through Novaseq and I'm not exactly sure how to trim my current RNASEQ data.

I've attempted to run the following script using trimmomatic

```
# go to working folder
cd /project/CHR_RNASEQ/reads

# run trimmomatic
trimmomatic SE -threads 8 -phred33 SLIDINGWINDOW:4:20 ILLUMINACLIP:$TRIM/adapters/TruSeq3-SE-2.fa:2:30:10 MINLEN:36 1_R1.fastq.gz 1_R1_trimmed.fastq.gz

```

It did log a job in the HPC, but I had no output at the end. I didn't receive an error either. This is unfortunately just some code I had lifted and trialed and errored until something finally ran, but I'm very sure I'm very wrong and any help with understanding this would be really appreciated.

Cheers

# Answer
> 1 votes
Depending on the cluster management tool, you might have received e-mails when the "job" begins and ends. If so, you can check the "Exit status" of the job. For example, in the case of our HPC the relevant lines from the "job completion e-mail" would be:

```
Execution terminated
Exit_status=0

```

BTW, at first I was confused with the phrase "Execution terminated", however, in the scheduler used in our system, it means "the run has been successfully or unsuccessfully completed". `Exit_status=0` int the next line means a successful run and for the unsuccessful runs there is a list with the error codes.

Secondly, if your run was successful, you would see new fastq files containing "\_trimmed" in the filename (reading this from your code). These fastqs should be more or less similar in size to the original depending on how aggressive the trimming settings are.

If you can provide more information, i.e. the type of cluster management tool used, the log file, ..., you can get more help and precise answers from the community.

---
Tags: rna-seq, data-preprocessing, trimming
---


# Question
Title: choosing a good representative genome subset
I'm trying to build a genomic database for DNA alignments. I started with NCBI accessions, but the data is very multiplicative, so I want to use subset of \[max\] N different strains for each specie. my question is: I can i obtain the minimal N that the samples will be diverse enough and will cover the specie good?

# Answer
Check out the progenomes database: http://progenomes.embl.de/
> 0 votes

---
Tags: sequence-alignment, clustering, subset
---


# Question
Title: Sequence alignment of extended IUPAC sequences
I want to align a short sequence to a reference, and both are encoded in the extended IUPAC notation (i.e. 16 letters).

What alignment software can I use for that?

For example, in Bowtie *"all non-A/C/G/T characters in reads become Ns",* so it ignores which extended letters are matching.

# Answer
> 1 votes
You can try using Biostrings, depends on how large is your reference, if it is small, you can just pairwiseAlignment with the nucleotideSubstitutionMatrix which supposedly :

```
library(Biostrings)
data(phiX174Phage)

x = DNAString("TCCAWRSCGCAGAAGT")

lapply(phiX174Phage,function(y)pairwiseAlignment(x,y,
substitutionMatrix = nucleotideSubstitutionMatrix()))

$Genbank
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$RF70s
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$SS78
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$Bull
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$G97
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$NEB03
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

```

If it is much larger, you can convert the sequence of interest to a pwm and find it in the reference fasta, like shown in this example

---
Tags: sequence-alignment
---


title: Finding gene length using ensembl ID

text:
# Question
Title: Finding gene length using ensembl ID
I want to find the length of a list of genes, of Homo sapiens, that is reported in the GEO database. I have gathered the ensembl id's of those genes.

I understand this information can be parsed from the start and end position given in GTF file available in ensembl.org.

I would like to ask for advice on how to programmatically download the gtf files with ensembl id's.

On the ensembl webpage a couple of options are listed, API and ftp access. I would like to know if sample scripts are available.

Update: I found a script here, it shows how to parse the data from gtf file. Could someone help in understanding how the code can be modified to obtain gene length from start and end position?

Edit 1: By start and end position, I refer to the start and end values given in the "fields" section here and the file from which I want to parse the values is the second last file available in this link.
For instance, the end - start value of the gene with ensembl id ENSG00000131482, gene :G6PC

Edit 2: I get the gene names and the expression values from the SOFT files documented in GEO, link

The following is a sample output from SOFT file. It gives me the probe\_set(ID\_REF) , Gene name and expression values from samples (GSMXXXXXX). Here we see fro G6PC, values are reported from two probe\_sets. Can we get the transcript id's from the probe\_sets?

Obtaining the transcript id will help me in finding the length of the specific transcript.

```
> Table(gds)[2449,0:5]
           ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
2449 1555612_s_at       G6PC   4.39472   4.52235   4.26869
> Table(gds)[16399,0:5]
         ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
16399 206952_at       G6PC   2.76693   3.45107    2.9815

```

# Answer
If you really just want the start and end positions of all genes in a list of genes from Homo\_sapiens.GRCh38.93.gtf.gz, you could do this in `awk`. For example, if your list of genes is:

```
$ cat genes
ENSG00000131482
ENSG00000223972
ENSG00000228794
ENSG00000187961

```

You can do:

```
$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$4,$5
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 11869 14409
ENSG00000228794 825138 859446
ENSG00000187961 960587 965715
ENSG00000131482 42900797 42913369

```

### Explanation:

* `if(NR==FNR){ wantedGenes[$1]++ }` : `NR` is the input line number and `FNR` is the current line number of *this* input file. The two are only equal while reading the first file. So, if this is the first file, save the first field of each line (the gene name) in the array `wantedGenes`.
* `if($3=="gene"){` : if this line is a 'gene' entry.
* `gsub(/";*/,"",$10);` : the gene names in the GTF file are quoted and have an `;` after the quote. For example: `"ENSG00000187961";`. Since we don't have this in the list, we need to remove it so we can match the gene name. This `gsub` command will remove all `"` followed by 0 or more `;` (the `;*` means 0 or more occurrences of `;`) from the 10th field of each line. So, it will replace `"ENSG00000187961";` with `ENSG00000187961`.
* `if($10 in wantedGenes){print $10,$4,$5}` : if the 10th field (the modified gene name from above) is in the `wantedGenes` array, then print the 10th field (the gene name), the 4th (start) and the 5th field (end).

You can easily modify this to print the length instead of the start and end positions:

```
$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$5-$4
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 2540
ENSG00000228794 34308
ENSG00000187961 5128
ENSG00000131482 12572

```

As requested in chat, here's a python version of the same thing:

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

wantedGenes = []
with open(sys.argv[1], mode='r') as f:
    wantedGenes = f.readlines()

wantedGenes = [ x.rstrip() for x in wantedGenes]

print('# Gene\tstart\tend\tlength')
with open(sys.argv[2], mode='r') as f:
    for line in f:
        if not line.startswith('#'):
            fields = line.split('\t')
            annotations = fields[8].split(' ')
            geneName = re.sub('[";]', '', annotations[1])
            if geneName in wantedGenes and fields[2] == 'gene':
                print('%s\t%s\t%s\t%d' %
                      (geneName, fields[3],fields[4],
                       int( int(fields[4]) - int(fields[3]) )
                      ))

```

Save the file as `foo.py` and run with the file names as arguments:

```
python foo.py geneList.txt Homo_sapiens.GRCh38.93.gtf

```
> 3 votes

# Answer
You could also do this in R with the GenomicFeatures library.

```
library(GenomicFeatures)

## make TxDb from GTF file
txdb <- makeTxDbFromGFF('Homo_sapiens.GRCh38.93.gtf')

## get gene information
all.genes <- genes(txdb)

## import your list of gene names
my.genes <- c('ENSG00000141510','ENSG00000184571','ENSG00000011007')

## get the length of each of those genes
my.genes.lengths <- width(all.genes[my.genes])
## put the names back on the lengths
names(my.genes.lengths) <- my.genes

## print lengths
print(my.genes.lengths)

ENSG00000141510 ENSG00000184571 ENSG00000011007
      25772           55687           18905

```
> 3 votes

# Answer
You can do this directly on BioMart site from EBI as well. Using a list of ENSG ids as input, you can get the actual length of each individual transcript, as well as the length of the genomic region coding said transcripts. These are the attributes I used:

```
Gene stable ID
Transcript stable ID
Transcript length (including UTRs and CDS)
Gene start (bp)
Gene end (bp)

```

The output is a table with the relevant data:

```
Gene ID Transcript ID   Transcript length (including UTRs and CDS)  Gene start (bp) Gene end (bp)
ENSG00000146648 ENST00000455089 3844    55019021    55211628
ENSG00000146648 ENST00000342916 2239    55019021    55211628

```

You can change the build as well if you don't want to use GRCh38.
> 0 votes

# Answer
If you are an R user, then you can easily do so with `biomaRt`

```
library("biomaRt")
library("dplyr")

genes <- read.delim("genes.file.txt", header=T)

head(genes)
          gene_id
1 ENSG00000000003
2 ENSG00000000005
3 ENSG00000000419
4 ENSG00000000457
5 ENSG00000000460
6 ENSG00000000938

## Build a biomart query
# In the example below, I use the human gene annotation from Ensembl release 82 located on "sep2015.archive.ensembl.org", More about the ensembl_build can be found on "http://www.ensembl.org/info/website/archives/index.html"

dataset="hsapiens_gene_ensembl"
mart = biomaRt::useMart("ENSEMBL_MART_ENSEMBL", dataset = dataset,
host = paste0("sep2015", ".archive.ensembl.org"), path = "/biomart/martservice", archive = FALSE)

annotations <- biomaRt::getBM(mart = mart, attributes=c("ensembl_gene_id", "external_gene_name", "start_position", "end_position"))
annotations <- dplyr::transmute(annotations, ensembl_gene_id, external_gene_name, gene_length = end_position - start_position)

# Filter and re-order gene.annotations to match the order in your input genes list
final.genes <- annotations %>% dplyr::filter(ensembl_gene_id %in% genes$gene_id)
final.genes <- final.genes[order(match(final.genes$ensembl_gene_id, genes$gene_id)),]; rownames(final.genes) <-NULL

head(final.genes)
  ensembl_gene_id external_gene_name gene_length
1 ENSG00000000003             TSPAN6       12882
2 ENSG00000000005               TNMD       15083
3 ENSG00000000419               DPM1       23688
4 ENSG00000000457              SCYL3       44636
5 ENSG00000000460           C1orf112      192073
6 ENSG00000000938                FGR       23213

```
> 0 votes

---
Tags: gene, ensembl
---


# Question
Title: Why index a sequence file with minimap2?
I have been using minimap2 for a few month now and I just found a script that follows these steps:

```
minimap2 [-x preset] -d target.mmi target.fa
minimap2 -a [-x preset] target.mmi query.fa > output.sam

```

I have never seen this `.mmi` index before. I always used minimap2 directly like this:

```
minimap2 -a [-x preset] target.fa query.fastq > output.sam

```

So my question is why should I index my reference genome before mapping? Is it only for speed up? Will it up speed my code if I do it? Or it is already done without specifying it?

# Answer
> 6 votes
If you don't provide an index, minimap2 will just generate the index on the fly (so, yes, it's done without specifying it). As such the mapping step will take longer. However, index generation is supposed to be very fast in minimap2 (a few minutes for the human genome), so it probably won't matter much.

Doing it separately will therefore also not speed up your code unless you end up using the pre-generated index multiple times.

---
Tags: minimap2, genome-index
---


# Question
Title: Optimum parameters for genome assembly via canu
I have been trying to assemble my genome via canu. Canu gave me about 200 contigs using default parameters. My genomesize is ~9Mbp. Can you recommend some parameters I can tweak to decrease the contigs number?

# Answer
The canu parameter reference guide is always a good place to start.

It would help to know a little bit more about what your assembly looks like. E.g.

* what is your n50, n90...?
* What is your dataset coverage?
* Have you run any tools like quast to evaluate quality?
* What is your input data type (probably PacBio)? Read length distribution?
* What is your expected chromosome/replicon number?
* Have you compared the assembly to any other data, e.g. genomes of close relatives? (quast can do this)

Basically, **why do you think that 200 contigs is too many contigs**? Sometimes you just don't have enough coverage to get more contiguous/"nice" assemblies. That's not to say it's not possible, of course, just that we don't have enough information to evaluate.

It's possible that some of those contigs are small and artefactual, so you might get rid of them with just length/coverage filtering. You can do some of this in canu itself (see "Output Filtering" section).
> 1 votes

---
Tags: assembly, canu, microbial-genomics
---


title: Extracting part of a csv by a list

text:
# Question
Title: Extracting part of a csv by a list
I have a .csv file like this

```
pair,S17A,S2..APOBEC.,S1..age.,S18.like,S3..BRCA.,S17B
LP2000104-DNA_A01_vs_LP2000101-DNA_A01,0.448027080410302,0,0,0,0,0.410602648202438
LP2000325-DNA_A01_vs_LP2000317-DNA_A01,0.343096214261112,0,0.111190317814775,0.171937435363758,0.187395051556141,0.156209735124657
LP2000327-DNA_A01_vs_LP2000319-DNA_A01,0.436249785689912,0,0.066118729867486,0.077392575088635,0.116384601626772,0.263251836199875

```

I have a list of samples which are a part of first column for example

```
LP6005334-DNA_H01
LP6008334-DNA_A03
LP6008334-DNA_A04
LP6008334-DNA_B02
LP6008334-DNA_C02
LP6008334-DNA_D02
LP6008336-DNA_H01

```

I want to extract .csv for this list of samples

I can not think of intersect function as these samples are a part of sample name not full name

Can you help?

# Answer
> 1 votes
You can use `strsplit`.

```
ciao <- "LP2000104-DNA_A01_vs_LP2000101-DNA_A01"
sampleIneed <- "LP6008334-DNA_A03"
sampleIneed %in% strsplit(ciao,"_vs_")[[1]]
[1] FALSE

"LP2000104-DNA_A01" %in% strsplit(ciao,"_vs_")[[1]]
[1] TRUE

```

You can loop (e.g. `sapply`or `lapply`) through the element of the first column to extract only the rows you need.

You can also simply use `grep` on the first column

```
ciao <_ your_csv
grep("LP2000104-DNA_A01", ciao$pair)

```

---
Tags: r
---


# Question
Title: Merge multiple PDB files into a single file
I have several pdb files and each of them corresponds to a domain on a protein. I want to merge these files into one single pdb for the whole protein, and the approach I'm looking for should allow me to order these domains according to the order by which they're located on the protein (for example, merging domain A, B, and C by the order CBA). It would be preferable if I can do the merging with VMD but any suggestion would be appreciated.

# Answer
> 1 votes
PyMOL is a rather flexible tool for doing this.

The command `create combo, domA or domB or domC` in the PyMOL app will create a single object. However, they will need to be the same chain and segment (`segi`). Also you will need to fix the numbering to match what you expect them to be. To fix this you need to do:

```
 alter domA, chain='A'
 alter domA, resv+=0
 alter domB, chain='A'
 alter domA, resv+=100
 alter domC, chain='A'
 alter domA, resv+=200
 alter all, segi=''
 sort
 create combo, domA or domB or domC
 save combo.pdb, combo

```

PyMOL can be installed as a Python 3 module, which works in the same way but with `pymol.cmd` and is a lot more replicable. But for starting out the App is fine.

```
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.alter('domA','chain="A"')
     etc.

```

Now. This assumed that the coordinates are where we want them to be. Which is highly unlikely and that there are zero overlaps.

* If there are overlaps, aligning the overlapping residues (`align domB and resi 80-100, domA and resi 80-100`) can be done. Just remember that the first element is the mobile object, the second the target, hence why I wrote it backwards.
* If do not have overlaps and want the domains docked together you can use the DEMO tool by the I-Tasser people.
* If you want your domains as "balls on a string" aligned on the x-axis for a nice picture and you do not care if it may not be biological, you may want to check out a Github repo of mine for doing just so (it offsets, aligns if possible and spaces them out, but does not add linkers).

If you want to add linkers, you can use submit your sequence to Phyre or I-Tasser with your missing density model as a reference, or you can add them with rosetta remodel (tricky) or make them up in Coot or PyMOL (bad idea, but if you are in a hurry... here is a tutorial of mine).

---
Tags: pdb
---


# Question
Title: a package with randomization/permutation tests for genomic analyses
please would you let me know, is there a package in R/python/C that has implemented **randomization/permutation tests for joint genomic analyses**

(an example of joint genomic analyses -- when jointly considering both GENE EXPRESSION and PROTEIN BINDING along the DNA).

The context of my question is the following :

let's consider **1000 UP\_REGULATED genes with increased PROTEIN X and with a HISTONE MARK Y :**

in order to show that PROTEIN X is related to HISTONE MARK Y for 1000 UP-REGULATED genes, what "controls" would you use for the comparison :

-- 1000 RANDOM GENES (and multiple randomization tests)

-- 1000 UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y

-- 1000 UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y

-- 1000 UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y

-- 1000 NOT-UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y

-- 1000 NOT-UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y

-- 1000 NOT-UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y

-- anything else ?

any suggestions are very welcome. thanks a lot,

bogdan

ps : 've posted a similar question on BioC website : https://support.bioconductor.org/p/128828/

# Answer
I would simply look at co-variance across the entire dataset and examine values \>0.8 in the first instance rather than select 1000 genes.

With your approach sampling 1000 genes can be criticised for selective sampling and the way to get around this is to perform replicates, with (bootstraps) or without replacement. You could produce e.g. 1000 replicates of 1000 genes and perform co-variance on each replicate to examine Y = X.

Python has a specific permutation/combinations module call itertools, which is different to random resampling and doesn't need replicates

```
from itertools import combinations
X = [1,2,3]
Y = [1,2,3]
print (list(combinations(iterable, 1000) + X + Y))

```

---

The standard way is to sample is,

```
from random import sample
mygenes = [1, 2, 3, 4, 5]

print(sample(mygenes,3))

```

**Output**

```
[2, 4, 5]

```

There is nuance and further funcationality ... I cannot remember is "sample" is with or without replacement (its quite a big deal)
> 1 votes

---
Tags: phylogenetics, statistics, phylogeny
---


# Question
Title: how to generate structure file from blast results of an input sequence
I have tried to incorporate blast in my application using biopython.i could extract the blast results but now i want to extract the structure file corresponding matches. My code is :

```
from Bio.SeqIO import parse
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

file = open("prats.fasta")

records = parse(file,"fasta")
for record in records:
    print("Id: %s" % record.id)
    print("Name: %s" % record.name)
    print("Description: %s" % record.description)
    print("Annotations: %s" % record.annotations)
    print("Sequence Data: %s" % record.seq)
    print("Sequence Alphabet: %s" % record.seq.alphabet)

from Bio.Blast import NCBIWWW
sequence_data = open("prats.fasta").read()
sequence_data

result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
result_handle

with open('results.xml', 'w') as save_file:
    blast_results = result_handle.read()
    print(blast_results)
    save_file.write(blast_results)

from Bio.Blast import NCBIXML
E_VALUE_THRESH = 1e-20
for record in NCBIXML.parse(open("results.xml")):
    if record.alignments:
        print("\n")
        print("query: %s" % record.query[:100])
        for align in record.alignments:
            for hsp in align.hsps:
                if hsp.expect < E_VALUE_THRESH:
                    print("match: %s " % align.title[:100])

```

from here i want to generate the xml files and then fetch structure of each nucleotides results

# Answer
## If you the DNA/RNA crystal structure:

Simply change the `nt` in `NCBIWWW.qblast("blastn", "nt", sequence_data)` to `pdbnt` to get the nucleotide sequence in a chain. The PDB code will be the accession ID with an underscore plus chain letter. Do note that the PDB structure may be offset.

## If you want the crystal structures of the encoded protein:

`blastx` (translation) on the `pdbaa` database. As above the PDB code will be as above.

## If you want a model of the encoded protein:

Check Swissmodel for a pre-existing close threaded model, else check out I-Tasser (best), Phyre (fast).

## If you want a model of the DNA/RNA:

You are running blast without special settings —its not a short motif—, so these things may be long, so do be weary of any result. There are a few RNA folding prediction tools online (e.g. RNAFold).
> 1 votes

---
Tags: python, blast, pdb, nucleotide-models
---


# Question
Title: How to calculate fold change in gene expression from RSEM (RNAseq)
I have gene expression data from RNAseq, specifically: log2(x+1) transformed RSEM normalized counts.

How can I convert these data from multiple samples to determine fold change in gene expression?

# Answer
You should use a proper statistical framework for RNA-seq dfferential analysis (which includes FC calculation). Standard tools for this are (among others) `edgeR` or `DESeq2`. You could use tximport to import RSEM outputs into R and then use its output for e.g. `DESeq2`. The linked manual provides example code for this. Please read the manuals of `edgeR` or `DESeq2` towards the respective downstream analysis. Both are extensive and compreshensive.
> 3 votes

---
Tags: rna-seq, gene-expression, rsem, fold-change
---


# Question
Title: HTSeq-count is returning 0 for every gene, instead of expression value
(Ported from stack overflow)

I'm trying to summarize gene count using htseq-count; and it's returning 0 counts at every gene. I'm not sure what I'm doing wrong; I think it has to do with the gene flag I'm using. I've tried using the GTF for Arabidopsis TAIR 10: which I got from: ftp://ftp.ensemblgenomes.org/pub/release-46/plants/gtf/arabidopsis\_thaliana/Arabidopsis\_thaliana.TAIR10.46.gtf.gz as well as the GFF file available on Araport.

I'm using various BAM files I know to have gene expression values.

I've tried both locally and on galaxy The local commands I've tried so far are:

```
python3 -m HTSeq.scripts.count --idattr=gene_id #bam #gtf > summary.txt
python3 -m HTSeq.scripts.count -i gene_id #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i locus_tag #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i ID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i SeqID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i Parent #bam #gff > summary.txt

```

I've done the same on galaxy; inputing the -i value in.

Any help as to where I'm going wrong would be great.

Thanks!

# Answer
Maybe the contig names in your BAM files don't match those in the GFF file? Ensembl uses contig names like `1`, some other reference sources use `chr1`.

Generally, for ensembl GFF files, HTSeq's default value for idattr should be fine though.
> 1 votes

---
Tags: sam, gene, gene-expression, gff
---


# Question
Title: What is the aim of insertion codes in the pdb file format?
In PDB file we always see insertion codes (ic) in a few proteins. Could anybody please tell me why it has been added in structure and what is the aim of that ???

# Answer
When I wondered it myself years ago and I found a good explanation in a mailing list archive. I think it was this one, by Frances C. Bernstein:

> The problem with numbering started when people wanted to compare the 'same' proteins from different species. They found that there were the following possibilities that gave rise to differences:
>
> 1. More or fewer residues at either end.
> 2. Extra residues at various places within the chain.
> 3. Fewer residues at various places within the chain.
> 4. Different amino acids at the same place.
>
> Now imagine that residue PHE 195 is very important for the activity of the protein in species A. But in species B it is residue PHE 197 and in species C it is PHE 212, because species B and C are not the same length as species A.
>
> Because people felt it was important to preserve the amino acid numbering for 'important' residues and to be able to readily discuss and compare the structures from different species, various people decided to try to number the proteins from species B, C, etc. to match the numbering used for species A. In doing this, one must have gaps (missing numbers) where a sequence is shorter. But what should one do when a sequence is longer? This is the case where it is necessary to insert extra numbers and this is done by using insertion codes.
>
> Thus the insertion code is an integral part of the residue number and it is improper to ignore that field when using a PDB entry. You must also allow for 'missing' numbers when using a PDB entry.

In the PDB format the insertion code is in column 27, following the sequence number (columns 23-26):

```
                      nnnnI
ATOM  11918  CZ  PHE D 100      -6.852  76.356 -23.289  1.00107.94           C
ATOM  11919  N   ARG D 100A     -9.676  74.726 -19.958  1.00105.71           N
...
ATOM  11970  CE  MET D 100H     -8.264  83.348 -19.494  1.00107.93           C
ATOM  11971  N   ASP D 101     -11.329  81.237 -14.804  1.00107.41           N

```

Together, columns 23-27 make a sequence ID, which is a mostly-numeric ID of the residue position.

Note: It's not a full ID of the residue. When you have a point mutation (a.k.a. microheterogeneity) you have 2+ residues with partial occupancy at the same position. To fully identify a residue both the sequence ID and residue name are needed.

As of 2020, the PDB (organization) uses PDBx/mmCIF as the primary format. The PDB format is only a legacy one. The presentation of the insertion codes (and many other things) in the mmCIF format are rather confusing, so it's worth explaining.

The mmCIF format introduced double identifiers for chains, residues and atoms: author's ID and "label" ID. Both look the same and are easy to confuse. For example: two chains named B and A by the author are renamed as, respectively, A and B by the wwPDB; both names are used in parallel: one chain is B(A), the other A(B). Poor design choice at every step. The same with the sequence ID. Originally, we had two tags:

* `_atom_site.auth_seq_id` -- author's sequence ID, possibly with an insertion code (100A),
* `_atom_site.label_seq_id` -- numeric-only sequence ID, numbered with respect to the sequence of this molecule (SEQRES).

Then, as I heard, numeric insertion codes were spotted in the wild. (In a database big enough any mistake that is not automatically checked for will happen). This made a problem with translation from the PDB format to mmCIF. The sequence number 100 with IC 1 would be interpreted as sequence number 1001. The numeric ICs must have been remediated in the meantime, because they absent as of 2020. But at that time (10-20 years ago?), instead of correcting the entries, the decision was made to move the insertion code to a separate column (pdbx\_PDB\_ins\_code). The documentation was never updated, so the auth\_seq\_id is still documented to store the sequence ID, not just the sequence number.

Here is an example from 3B9F:

```
                             !!  I                                      nn
ATOM   770  C CD2 . LEU B 2  46  ? -23.416 21.041  -27.401 1.00 19.94 ? 60   LEU H CD2 1
ATOM   771  N N   . TYR B 2  47  A -21.191 16.717  -28.119 1.00 18.73 ? 60   TYR H N   1

```

The sequence IDs here are 60 and 60A. Since the meaning of the insertion code is not widely known, it landed next to label\_seq\_id which has nothing to do with it. So again: the insertion code is just a part of the author's sequence ID (originally, in auth\_seq\_id), and the letter A above (10th field) is simply an extension to the number 60 (17th field).
> 7 votes

---
Tags: pdb
---


title: using blastn to search [error]

text:
# Question
Title: using blastn to search [error]
so I'm fairly new to using linux and blast - have encountered an error in my work.

I have created a fasta database with the results in this directory path: `~/class/hw/db`

I'm currently trying to run a search against my created database against query sequences in a file called `seqs.fas`.

My code is as follows: `blastn -db $db -query seqs.fas -out q1results.txt`

However I keep getting this error: `Error: Too many positional arguments (1), the offending value: seqs.fas`

`Error: (CArgException::eSynopsis) Too many positional arguments (1), the offending value: seqs.fas`

I'm not sure what I'm doing wrong here, and would appreciate any advice to fix this, thanks!

# Answer
I remember having this problem in the past, it was encoding related. Rewriting the `-`s manually (the dash sign) should fix the issue, please see this.

> Problem is a copy & past from the BLAST documentation, all dashes ('–') need to be replaced by hand at the command line.
> 2 votes

---
Tags: blastn
---


title: Convert PAF format to SAM/BAM format

text:
# Question
Title: Convert PAF format to SAM/BAM format
I have a bunch of PAF files resulting from the alignments of fastq files on a reference genome with minimap2. I would like to convert them into SAM/BAM format so I can use `samtools mpileup` to recover the read coverage at filtered positions

So far I've looked at `paftools.js` and found `sam2paf` but no `paf2sam`

I am open to suggestions but I'd like to avoid mapping the reads again.

**Edit**

I have opened an issue on minimap2's github

Answer from Heng Li :

> You need to rerun alignment. No tools exist to convert paf to sam.
>
> ```
> why there is no paf2sam in paftools.js
>
> ```
>
> It is technically possible, but someone needs to write the code.

# Answer
> 2 votes
The paf format does not have the CIGAR string and base quality information that are typically present in a BAM/SAM file. While it is in principle possible to write a paf2sam converter, it's going to be very difficult to regenerate the CIGAR string with all the indel information (and maybe it's probably impossible to do so) that you'll need in a BAM/SAM file. If your motivation is to get coverage information, it might make more sense to convert the paf file into a bed file, before using a tool like bedtools to convert it into a coverage track.

---
Tags: file-formats, format-conversion, sam, minimap2
---


# Question
Title: Where can I find summary data for how common certain mutation *types* are?
I'd like to know how common certain mutation *types* are in public data sets like the 1000 Genomes, ExAC, and ESP6500. Specifically, I'd like to know the distribution of stop-gains, stop-losses, frameshift, and other mutation types. For example, what is the median number of stop-gains observed across individuals in the ExAC population? I'd like to have the full distribution, but median values would work.

To be clear, I am *not* looking for allele frequencies.

I didn't find this kind of summary data when I looked at the respective websites. Is there a good source for these summary data?

# Answer
This dbsnp github repo seems to offer that information.

https://github.com/ncbi/dbsnp/tree/master/Build%20Announcements/151
> 1 votes

---
Tags: public-databases
---


# Question
Title: How to construct the phylogenetic tree of the whole genus?
I am confused about how to construct the phylogenetic tree of the streptomyces genus. The taxonomy of the streptomyces genus is mostly based on the MLST method, but now I need to use the gene data from NCBI to construct the phylogenetic tree of the streptomyces genus. Should I use the whole gene of the MLST's housekeeping gene to construct the tree? Or I should just use the loci to build the tree? Or the other method? I am a junior in university, have little knowledge of the filed., but my tutor asked me to do that. plus: Sorry for my awkward expression~

# Answer
The MLST genes are not used to contruct a tree, because each pattern of alleles determines a predetermined genotype, e.g. Staphlococcus aureus ST 22 is the MLST type concerning the MRSA epidemic. You could build a tree, but its only purpose is to visualise the ST, so its a bit pointless but okay at least you will not mistype something.

Performing a whole genome tree on Strep for the genus is complicated but doable and there are missing genes and stuff and if you just want typing then its not needed.

---

So if typing is your goal I would submit your MLST loci to the Steptomyces MLST database,

https://pubmlst.org/bigsdb?db=pubmlst\_streptomyces\_seqdef
> 1 votes

---
Tags: phylogenetics, phylogeny
---