# Question
Title: Confusing result from Sanger Imputation Service (Eagle v2.4 for phasing, PBWT v3.1 for imputation)
I ran my 23andMe data through the Sanger Imputation Service that uses Eagle v2.4 for phasing, PBWT v3.1 for imputation.

However, some aspects of the results are very confusing to me. Perhaps see #3 (the most confusing) first.

**Point 1**

1. There are 62,623 rsIDs for chromosome 3 in my 23andMe file and 2,821,894 SNPs (actually, unique rsIDs) after imputation. However, only 60,144 of the original 62,623 rsIDs are in that ~2.8 million. Where did the other 2,479 go? Here are some examples of missing rsIDs: rs1516332, rs12374025, rs1170695, rs2731343 and rs63749864.

**EDIT Point 1**

Looking at `rs1516332`, the position is 3:105793, which does have a result in the file, but it doesn't have an rsID:

```
3    105793    .    C    T    .    PASS    RefPanelAF=0.17878;AN=2;AC=0;INFO=1    GT:ADS:DS:GP    0|0:0,0:0:1,0,0

```

Looking at `rs63749864` (3:13564048) it's missing by rsID and position.

**EDIT Point 1 (again)**

As described here: https://github.com/richarddurbin/pbwt/issues/51

> "... of the 2,479 'missing' SNPs, 1,806 can be 'found' \[by position\] leaving 673 'anchors' \[or rsIDs called by 23andMe\] missing from chromosome 3."

**Point 2**

2. Some imputed rsIDs appear more than once in the results. How should I interpret that? e.g.

```
3   64071   rs116059446 G   T   .   PASS    RefPanelAF=0.000369572;AN=2;AC=0;INFO=1 GT:ADS:DS:GP    0|0:0,0:0:1,0,0
3   64071   rs116059446 G   A   .   PASS    RefPanelAF=9.2393e-05;AN=2;AC=0;INFO=1  GT:ADS:DS:GP    0|0:0,0:0:1,0,0

```

**EDIT Point 2** From the answer below, this is a variant with more than two alleles: rs116059446

Seems the format of this VCF file only allows one alt per line, hence the two alts are split over two lines.

See my question here: Is it valid VCF not to 'squash' positions with more than one ALT allele?

And the answer from Durbo here: https://github.com/richarddurbin/pbwt/issues/52#issuecomment-1062951105

In the above case, the genotype is easy to interpret (REF/REF in both lines).

However, I'm now finding examples that I just can't interpret:

```
3   255395  rs331869    C   G   .   PASS    RefPanelAF=0.385864;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    1|0:1,0:1:0,1,0
3   255395  rs331869    C   T   .   PASS    RefPanelAF=0.980998;AN=2;AC=2;INFO=1    GT:ADS:DS:GP    1|1:1,1:2:0,0,1

```

*Is my genotype CG or TT at this position?*

**Point 3**

3. Sometimes a rsID appears twice but at different locations:

```
3   4942430 rs71634747  G   A   .   PASS    RefPanelAF=0.445673;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    1|0:1,0:1:0,1,0
3   4942432 rs71634747  C   G   .   PASS    RefPanelAF=0.248907;AN=2;AC=0;INFO=1    GT:ADS:DS:GP    0|0:0,0:0:1,0,0

```

Unless I'm really confused, this just seems wrong. If it wasn't for this weirdness, I'd seek concrete answers for #1 and #2, but it could be that they are both just a symptom of this potential error?

# Answer
> 3 votes
1. Don't know. These could be retired or inconsistent SNPs (e.g. shifted locations, or otherwise found to be different from expected). It might be that the imputation suggests that the SNP should be something other than what is observed. I can't see anything obvious from the dbSNP information that I've looked at.
2. a) In this specific case for rs116059446, that's a SNP that has more than two variants. One way of normalising VCF files is to rearrange them so that each line only refer to a single polymorphism from the reference sequence, so if there are two non-reference variants, the variant will get two lines in the VCF file.

b) For rs331869, this is another SNP with more than two variants recorded in dbSNP. Given that you're digging through the weeds of a huge dataset, the way of thinking about this needs to change from "Is this common" to "Is this possible". Problems are much more likely to exist with SNPs that aren't biallelic, because it indicates that something funky is going on with genetics or testing. I'll propose three possibilities for what could be happening here, in what I consider to be descending order of likelihood:

1. you have one non-reference allele on one of the sister chromosome copies, and another non-reference allele on another copy
    2. you have a somatic mutation (i.e. a mutation that happened during your lifetime), and some of your cells have one variant, and others have a different variant.
    3. the probes are tagging different regions in your genome from what is expected. Probes should be location-unique, but that's impossible to guarantee for everyone, because genetic variation is a thing that happens.
3. In this specific case for rs71634747, it looks like that SNP has been retired, and has merged with rs724135. It's possible that this is a multi-nucleotide variant where variants are only observed in tandem. This is an intronic variant, so seeing a 3-nucleotide variation is a little unexpected, but given that you're picking out exceptions from a very large dataset, not completely unexpected.

I think the main take-home from this is that you shouldn't ever expect results to be a perfect representation of reality. It can be useful to look at raw data to help understand what's going on, but there's too much weird stuff going on with biology, chemistry and physics for anything to be 100% accurate.

If you did exactly the same test at a later date, some of your 23andMe results will be different. If you did imputation at a later date, with a different program version, or different imputation reference dataset, your results are likely to be different.

Accuracy can be substantially improved by using different technologies to attempt to discover the same thing, and only concentrating on consistently-reported differences. You could try an Infinium Omni2.5 SNPchip, for example, or whole-genome sequencing. Unfortunately, this will lead to another problem: what should be done with the inconsistent results? The answer to *that* is situation dependent, and is one of the hardest challenges of bioinformatics.

# Answer
> 0 votes
The answer is to relax and throw away the 0.001% of stuff that doesn't make sense, lower your standards and just get on with finishing the job you set out to do without actually caring about consistency, logic, meaning or formats that are actually defined or interpretable.

Nobody cares about what these SNPs mean, how to interpret them programmatically, why there are version issues in the output of the imputation server or why it throws away ground truth in the output, leaving the user to struggle with all these issues. I mean... if nobody cares, why should I?

I've written up as detailed a description of the problems as I can here: https://github.com/richarddurbin/pbwt/issues/52 and here: https://github.com/richarddurbin/pbwt/issues/51

and written a script to throw away all the junk here: https://github.com/Geromics/covcheck/blob/wip/report-v2/Research/debug\_imputation\_results.py

---
Tags: vcf, snp, imputation, pbwt
---


# Question
Title: Can you style network visualizations of directed graphs using NETWULF?
I want to draw directed graphs using networkx, followed by manual tuning of the visualization with the use of NETWOLF.

Here is the example of a network.

```
G =nx.DiGraph()
G.add_edges_from([(1, 2), (1, 3),(1, 4),(2, 3)])
network, config = nw.visualize(G,config={'zoom':3})

```

However, the graph shown in NETWOLF was undirected graph.

Do you know if you can style network visualizations of directed graphs using NETWULF?

# Answer
I'm the author of gravis, an open-source package for interactive graph visualization in Python. It allows to draw directed graphs from NetworkX and the resulting visualization can be modified by user interaction (e.g. node positions, edge curvature and various other visual properties).

Here's a minimal example based on your code:

```
import gravis as gv
import networkx as nx

g = nx.DiGraph()
g.add_edges_from([(1, 2), (1, 3), (1, 4), (2, 3)])
gv.d3(g)

```

Output: NetworkX DiGraph visualized with gravis inside a Jupyter notebook:
> 3 votes

---
Tags: python, networks, graphs
---


# Question
Title: How to make a network overlap of gene lists across 3 conditions like this image?
From 3 gene lists (lung, nasal, gut), how to make this gene overlap network plot?

This image is figure 4a from: https://www.nature.com/articles/s41591-020-01227-z

# Answer
I'm the author of gravis, an open source package for interactive graph visualization in Python.

* It can visualize directed and undirected graphs (=having edges with or without arrows) as well as multigraphs (=having more than one edge between a pair of nodes).
* It allows to control visual properties such as node sizes, colors, shapes or edge widths and color with data.
* It recognizes graph objects from various network analysis packages such as NetworkX, igraph, graph-tool and others.
* The output is based on web technologies (HTML/JS/CSS), so it can be viewed in a browser or Jupyter notebook. There is also support for exporting static images (JPG, PNG, SVG).
* It comes with many examples on the documentation site.

Here's a quick attempt to reproduce a part of the figure you've provided:

```
import gravis as gv
import networkx as nx

s1 = 40
c1 = 'white'

s2 = 20
c2 = '#ffd703'

s3 = 20
c3 = 'orange'

s4 = 20
c4 = 'red'

g = nx.Graph()
g.add_node('Lung', color=c1, size=s1, border_size=1)
g.add_node('Nasal', color=c1, size=s1, border_size=1)
g.add_node('Gut', color=c1, size=s1, border_size=1)

g.add_node('MDFI', color=c2, size=s2)
g.add_node('LOXL1', color=c2, size=s2)
g.add_node('CNBP', color=c2, size=s2)

g.add_node('C3', color=c3, size=s3)
g.add_node('CA2', color=c3, size=s3)
g.add_node('SEPP1', color=c3, size=s3)

g.add_node('TFF3', color=c4, size=s4)
g.add_node('PTMA', color=c4, size=s4)

g.add_edge('Lung', 'MDFI')
g.add_edge('Nasal', 'LOXL1')
g.add_edge('Gut', 'CNBP')

g.add_edge('Lung', 'C3')
g.add_edge('Nasal', 'C3')

g.add_edge('Lung', 'CA2')
g.add_edge('Gut', 'CA2')

g.add_edge('Nasal', 'SEPP1')
g.add_edge('Gut', 'SEPP1')

g.add_edge('Lung', 'TFF3')
g.add_edge('Nasal', 'TFF3')
g.add_edge('Gut', 'TFF3')

g.add_edge('Lung', 'PTMA')
g.add_edge('Nasal', 'PTMA')
g.add_edge('Gut', 'PTMA')

gv.d3(g)

```

Resulting output in a web browser:

I've moved the nodes manually a bit so that the layout resembles the provided image. Once everything looks as intended, a static image can be exported with a single button click on the right-hand side.

There are many other ways to produce static graph visualizations, e.g. programmatically with other Python libraries mentioned before, or interactively in standalone tools such as Gephi, Cytoscape or Tulip.
> 1 votes

# Answer
You could use OmniGraffle to make this figure.

Objects in OmniGraffle have so-called magnets to which you can attach lines, which can join two objects.

Once you have objects joined, you could group them and position them as desired.

As you drag an object, any attached lines will also change their position as needed. You can drag groupings of objects (e.g., conditions or other categories of objects), which does the same thing.

You can export an illustration to PDF or SVG format from OmniGraffle, which will let you do any further work in Adobe Illustrator or other vector illustration tools, or share the illustration with colleagues over the web, for instance.
> 0 votes

# Answer
This question is a bit difficult to fathom, but essentially I think the answer you are look for is called network theory. Its quite simple, but alas not fashionable. You'll find an introduction here and try and avoid 'directed graphs' because that gets complicated.

The alternative way is to use Google fusion tables, which incorporates network analysis, which is really easy to use and investigate. Its very interactive. You'll find more information here, https://support.google.com/fusiontables/answer/2566732?hl=en

Thats it and nothing more, undirected network theory is pretty simple. Its does not tell you very much (in my opinion), but it makes pretty diagrams.

Please note there is an analysis called 'Graph theory' and is not needed here. Whilst this has some similarity to 'network theory' it is distinct and the two shouldn't be confused.
> 0 votes

---
Tags: visualization
---


# Question
Title: Visualize network (graph) using "line" layout
I want to visualise a network in a way that the nodes are along a line (and the edges are not on that line). I have included a simple example below.

Which network visualisation programs (Cytoscape, Gephi, R igraph etc) do have this layout option?

# Answer
That plot is called an Arc diagram you can google that and use the package or program that best suits you.

Here are 2 examples of the use of Arc diagrams in R:
Example 1:

Example 2:
> 2 votes

# Answer
I'm the author of gravis, an open source package for interactive network visualization in Python. It supports many ways to influence visual properties with data, such as node size, color or location. It also allows to control edge curvature. Therefore it's possible to place all nodes along a line (or another form) and have edges with a chosen degree of curvature between them.

Edit: Python code that reproduces your image with NetworkX and gravis

```
import gravis as gv
import networkx as nx

g = nx.Graph()
g.add_node(1, x=0, y=0)
g.add_node(2, x=50, y=0)
g.add_node(3, x=100, y=0)
g.add_node(4, x=150, y=0)
g.add_edge(1, 2)
g.add_edge(1, 3)
g.add_edge(3, 4)
g.graph['node_color'] = 'blue'
g.graph['node_size'] = 15
gv.d3(g, edge_curvature=1.5, show_node_label=False)

```

Result inside a Jupyter notebook (optional, just for convenience):
> 2 votes

---
Tags: networks, visualization
---


# Question
Title: What is the basic difference between a protein and a ligand?
Is ligand a type of protein or vice versa?

# Answer
> 1 votes
A ligand is any molecule that can bind to a receptor. A receptor is always a protein. A ligand can be a protein, but also something else, be it a combination of protein with carbohydrates and-or fats, or just molecules such as cyclic AMP, any chemical compound, structures expressed on the cell surface, pretty much any structure you can think of.

---
Tags: proteins
---


# Question
Title: Can I use samtools addreplacerg to replace multiple RG entries at the same time?
I have a bam file that contains two @RG lines:

```
@RG     ID:Foo.Bar_L1.Foo   PU:Bar_L1.Foo SM:Foo        PL:Illumina     LB:Library.fa
@RG     ID:Foo.Bar_L2.Foo   PU:Bar_L2.Foo SM:Foo        PL:Illumina     LB:Library.fa

```

I wish to replace it with one @RG line:

```
@RG     ID:Foo SM:Foo        PL:Illumina     LB:Library.fa

```

What is the best way to do this?

samtools addreplacerg looks useful, but it looks like I can only affect one @RG per call, and I have no desire to write out the whole bam file twice

# Answer
> 2 votes
Not sure about samtools, but picard AddOrReplaceReadGroups can clear existing read groups and write a new one.

---
Tags: sam, samtools
---


# Question
Title: Heatmap of motifs present in protein sequences
I have a file containing multiple protein amino acid sequences. I am interested in visually presenting the location of amino acid motifs PQG, QQG and AQG in the sequence in the form of a heatmap in R. Every location of the motifs should be represented as a differently coloured barcode. Example sequence follows,

```
>cat archae.fasta
>A0A075FKY9
MRRITLPKGIEISGKSPFGTSSGVYEVVISISENKPSDASIASNSFNSLWKDNFHLRVSG
DHFSQVLGSGSNPIPDNVFDLGSVWIVIQDQFSPVHTSFQFNISPSNQPASKPEEKPRRV
ETSTVKRTITKPGRLGPPGGKGPRGPSGYPGIKGDKGGRGPTGDKGDKGDKGIPGPQGEK
GKTGQAGDKGDKGITGPPGEKGPRGSTGPPGDKGDKGVQGKIGDKGLTGTTGPVGDKGTQ
GPQGPPGGKGLTGIPGPQGDKGTKGPLGPVGERGPTGKPGESGLQGPQGIQGAQGERGHS
GPPGPQGEKGLIGDKGEIGERGTRGPPGPPGEKGSQGGMSDEGKRLIKELLELLASKNII
STEEQIKLTSYLY
>A0A075FMB0
MPGDKGTSGIRGVPGDKGDKGPQGPPGDKGLTGPVGVPGEKGPTGTSGVQGEKGIQGTPG
PLGEKGVTGPAGDKGPIGPPGPAGSKGITGPAGPLGDKGTIGPPGPLGDKGSKGPEGPTG
DKGPQGPQGPAGAKGLTGVPGPQGEKGEKGPLGPIGEKGTTGATGPPGDKGPQGPQGTQG
ERGTTGPSGEKGPQGPQGIQGPQGERGPTGPIGSIGETGPIGGQGPVGPAGPRGPPGPPG
EKGPSGGMSSEQKAIFKELLEILTTKEIISTEEQIKLMSYLY
>A0A075G2F8
MPKGIEISGKSPFGTSSGVYEVKISVSQSKPSNDSITNNSFNSLWKDNFHLRVSEGYFSE
ILGSDSNPIPDNVFDLGSVWIVIQDQFSPVHTSFEFKISQDSKSIPTPKEKPRRVDTSSV
KKVRTIPTRQGPTGSKGIPGERGYPGTKGEKGGKGPTGDKGDKGDKGVPGPQGEKGKTGS
TGDKGDKGITGSPGEKGPRGPTGPPGDKGDKGVQGKIGDKGLTGTTGPIGDKGPQGPQGP
PGDKGLTGIPGPQGDKGGRGTIGPSGERGPTGRPGESGLQGPQGIQGPQGERGHSGPPGP
QGEKGIIGDKGEIGERGPRGPPGPPGEKGSQGGMSDDGKILFKELLQLLASKNIISTEEQ
IKLTSYLY
>A0A075GLR7
MPKGIEISGKSPFGTSSGVYEVRISISQDKPSDDGITNNSFNSLWKDNFHLRVSEGHFSE
ILGTDSNPIPADVFGLGSIWIIIQDQFSPVHTSFELHISQDNKLLSVPQGKPKPVNATKV
QKVRAQPGRPGPQGEKGGKGPAGYPGTRGEKGGRGITGDKGDKGDKGIPGPQGEKGIKGP
AGDKGDKGITGPVGPSGPKGSTGPPGDKGDKGVQGKSGDKGLTGVTGPLGDKGTQGPQGP
PGPKGLTGVPGPQGDKGMRGSIGPAGERGPTGKPGDPGLPGPQGIQGPQGERGHTGPQGP
PGEKGLVGDKGEIGDRGPRGPPGPPGEKGSQGGMSEESKRLVKELLELLASKNIISTEEQ
IKLTSYLY

```

As an example, please see a two-colour heatmap. I need something similar but with three colours corresponding to each of the three motifs.

I looked at the heatmaply package but the biggest problem is that my matrix is not numeric. Is there way to make heatmap using text strings in R?

Motif is defined as any three character string that matches the general pattern either PQG, QQG and AQG in the uploaded sequence. y-axis is the fasta code appended after \> (A0A075FKY9 for example). x-axis would depend on how the data is converted into a matrix.

# Answer
> 0 votes
Firstly, the problem with your data is that it ain't aligned. Without that your motif positions are not relative from one protein to another. Basically you need to align it and there's no wriggle room.

Secondly, its easy you simply parse the motifs you want and their location within the string - after it is aligned. In Python you would use `finditer` within a comprehension. This will go into a dataframe and subject to a heatmap package, I'm not a 'heatmap' person but it would be easy to graph where the x axis is the position where each motif is found and a qualitative score given for each type of motif (1-3 here). I don't know the equivalent in `R`.

There is an issue whether an indel will remove a motif and you would need to check for that by running the same command against an unaligned protein and seeing if the frequency matches.

My advice is to ask a new question with some code. I dunno whether someone has coded this already, - its easy to do. Personally, I don't think there's sufficient requirement, but I could be wrong.

---
Tags: r, heatmap, motifs
---


# Question
Title: Assemble Anaeroplasma species genome from metagenomic PacBio data
I have a fasta file containing reads generated by PacBio HiFi whole genome sequencing of a feces sample from mouse.

I would like to use this dataset to generate an assembled circularized genome for an Anaeroplasma organism known to be in this sample.

I am familiar with Linux and R. I would very much appreciate your advice on the needed steps and the required tools or R packages.

# Answer
> 1 votes
Try hifiasm-meta (paper here). Whether it works depends on the abundance and the complexity. You may also try HiCanu.

# Answer
> 0 votes
I would suggest the following workflow, understanding you to have metagenomic HiFi sequencing that you **know** to include your species of interest at high abundance:

1. metagenomic assembly of HiFi with metaFlyE
2. Annotate your resulting long contigs (\>500Kbp) using GTDB-tk
3. Look for your Anaeroplasma among long contigs
4. Assess whether your genome is closed.
5. If the genome is closed, you are done.
6. If the genome is not closed, you need to close it. Assess its completeness using e.g. checkM.
7. Map your HiFi reads to the genome, see if you can close the remaining gap with e.g. SSpace-LongRead.

I'd suggest using this paper as an example for further workflow details.

---
Tags: assembly, linux, wgs, pacbio
---


# Question
Title: How can I download from NCBI all the ITS genes and the related taxonomy?
I would like to download all the ITS1 and ITS2 genes from NCBI in a fasta file. And, I'd like to download even the related taxonomy of each sequence.

Thanks,

Marco

# Answer
> 4 votes
I recommend to use **Biopython**:

```
from Bio import Entrez

Entrez.email = "your@somemail.com"
Entrez.api_key = "YOUR_API_KEY"

search_result = Entrez.read(Entrez.esearch(db="nucleotide", term="ITS1[Gene Name]"))
with Entrez.efetch(
    db="nucleotide", id=search_result["IdList"], rettype="fasta", retmode="text"
) as handle, open(filepath, "w") as file:
    file.write(handle.read())

```

# Answer
> 3 votes
For batch operations, it would be better to acquire a NCBI API key and combine GNU parallel and Entrez Direct . If you already have the accession numbers of ITS1 and ITS2 genes, put the accession numbers in a text file like "accession-number.txt". Otherwise search the ITS genes you want to download and then get the accession number:

```
esearch -db nuccore -query 'you query term[Gene Name]' |
  efetch -format acc >accession-number.txt

```

To download all the sequences in fasta format.:

```
parallel -a accession-number.txt -j8 'efetch -db nuccore -id {} -format fasta >{}.fa'

```

Here `-j8` is used to control the number of jobs run in parallel, it should be less than the limit of 10.

To get related taxonomy (for eg, taxonomy rank, scientific name, taxonomy id, genus and species) of each sequence:

```
body='BEGIN { OFS = "\t" } { print "{}", $0}'
paralle -a accession-number.txt -j8 "elink -db nuccore -id {} -target taxonomy -name nuccore_taxonomy |
  esummary |
  xtract -pattern DocumentSummary -element Rank ScientificName TaxId Genus Species |
  awk '$body'"

```

---
Tags: phylogenetics, taxonomy
---


# Question
Title: How can I get latest .vcf files with annotation data?
I'd like to perform the annotation (1000genomes, COSMIC etc.) of my variants using SnpSift and SnpEff, however so far all I get are vcf files for separate chromosomes: http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/

I'd like to annotate my variants and then filter out common SNP's. What in your opinion is best annotation source and where can I get the file?

# Answer
> 2 votes
Do you mean `Somatic Sequence Variants`?,

if so, Table 1 inside the following full-text link: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5707196/

And additional information from my search--\> Here's https://docs.google.com/document/d/1hVtgMYAtvDoUcigKx8OLDmongDRzDQAd/edit?usp=sharing&ouid=102408774767674945852&rtpof=true&sd=true

# Answer
> 1 votes
According to your answer I tried to get the files which might be useful regarding annotation:

---
Tags: vcf, sequence-annotation, ngs
---


title: gtf2bed only retains gene information

text:
# Question
Title: gtf2bed only retains gene information
I want to convert an hg38 GTF file to bed format using `gtf2bed`:

```
gtf2bed Homo_sapiens.GRCh38.106.gtf > Homo_sapiens.GRCh38.106.gtf.bed

```

It outputs a bed file but it only contains the coordinates for the whole genes and has lost the locations of the exons, UTRs, etc. While I am running `gtf2bed`, I get thousands of warning like this:

```
WARNING:root:None   21  33276625    33276753    +   three_prime_utr
WARNING:root:WARNING End is after gene end (33279918 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  33279752    33279918    +   three_prime_utr
WARNING:root:WARNING End is after gene end (33283212 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  33283094    33283212    +   three_prime_utr
WARNING:root:WARNING End is after gene end (37101343 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37100814    37101343    +   transcript
WARNING:root:WARNING End is after gene end (37101343 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37100814    37101343    +   exon
WARNING:root:WARNING End is after gene end (37237744 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37221419    37237744    +   transcript
WARNING:root:WARNING End is after gene end (37221597 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37221419    37221597    +   exon
WARNING:root:WARNING End is after gene end (37237744 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37237404    37237744    +   exon

```

How can I create a bed file that also contains the coordinates of the other genomic features?

# Answer
I cannot reproduce this. Be sure to use the latest version of BEDOPS and to use the correct command:

```

wget http://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.chr.gtf.gz

gtf2bed < <(zcat Homo_sapiens.GRCh38.106.chr.gtf.gz) > Homo_sapiens.GRCh38.106.chr.bed

$ head Homo_sapiens.GRCh38.106.chr.bed
1       1211339 1214153 ENSG00000186827 .       -       ensembl_havana  gene    .       gene_id "ENSG00000186827"; gene_version "11"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding";
1       1211339 1214153 ENSG00000186827 .       -       ensembl_havana  transcript      .       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1213982 1214153 ENSG00000186827 .       -       ensembl_havana  exon    .       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; exon_id "ENSE00001832731"; exon_version "2"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1213982 1214127 ENSG00000186827 .       -       ensembl_havana  CDS     0       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; protein_id "ENSP00000368538"; protein_version "3"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1214124 1214127 ENSG00000186827 .       -       ensembl_havana  start_codon     0       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";

```

Finishes without any of the mentioned warnings.
> 3 votes

---
Tags: bed, gtf, gffutils
---


# Question
Title: Convert Arlequin ("DNA" markers) to Genepop (PGDspider or otherwise)
I have Arlequin output files (\*.arp) from fastSimcoal2 that I'm trying to convert to genepop files (to read into adegenet). The Arlequin files are using the "DNA" marker (50 bp long), and each includes ~2,500 samples.

When trying to use the program PGDspider to convert these (either the GUI or the CLI; v2.1.1.5), things proceed smoothly until the 1,000th allele is reached, where the program generates an error:

```
INFO  07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:596)  -Allele "TGGCATCAGCTTCTTCTAAGCAGGCCCGCTTTGTCTAGCGACAGATCTCG" converted to "999"
INFO  07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:596)  -Allele "ATCTTGTGGTCATGGAGCACCCAGAGCCAAACTCCTAAAGAAGACCCCGC" converted to "1000"
ERROR 07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:611)  -data "1000" in ind "2_375" (population: Sample 2) is too long!!!
INFO  07-03 11:18:10,631 (GenepopWriter.java:parseAlignedSampleData:618)  -In GENEPOP, alleles cannot be coded with more than 3 digits.

```

Is there a way around this error? I've processed RADseq data with many 1,000s of alleles using adegenet before, so I know these types of files can be written into genpop--I just don't know how to specify this to PGDspider. There's a field in the .spid file that I might be able to use to fix this (see below), but I don't know what to populate that field with, and there aren't any suggestions from the (limited) documentation.

```
# Specify the locus/locus combination you want to write to the GENEPOP file:
GENEPOP_WRITER_LOCUS_COMBINATION_QUESTION=

```

Alternatively, if there are suggestions for running this conversion without using PGDspider, I'm all ears.

Thanks for your help!

# Answer
> 2 votes
Ultimately, I achieved this using the R package strataG. The steps taken (in R) are shown below.

```
arp.object <- arlequinRead(arp.path)
gtype.object <- arp2gtypes(arp.object)
genind.object <- gtypes2genind(gtype.object)

```

I found that fastSimcoal2 shows an odd behavior when simulating DNA markers, in which sequences are concatenated together (regardless of marker blocks specified; relevant post here). Running the fastSimcoal2 simulations from strataG addressed this bug (and is more straightforward/reproducible, anyways).

---
Tags: linux, phylogenetics, format-conversion
---


# Question
Title: Generate a single column from two different columns in the same format
I have multiple data frames in R in the following format:

`Dataframe1`:

```
Var1    Var2    cor p   fdr
Bacteroides_galacturonicus  Lactobacillus_rogosae   0.967374392060265   0   0
Bacteroides_galacturonicus  Lachnospira_pectinoschiza   0.526779994497013   0   0
Lactobacillus_rogosae   Lachnospira_pectinoschiza   0.29419393866076    0   0
Bacteroides_galacturonicus  Bacteroides_salyersiae  0.374581964109831   0.548   0.829267117726658
Lactobacillus_rogosae   Bacteroides_salyersiae  0.45736458243079    0.002   0.03622024291498
Lachnospira_pectinoschiza   Bacteroides_salyersiae  0.566636776534  0.21    0.56910941475827
Bacteroides_galacturonicus  Holdemanella_biformis   0.342443352 0.094   0.394447279549719
Lactobacillus_rogosae   Holdemanella_biformis   0.15592768883205    0.174   0.527328455284553

```

`Dataframe2`:

```
Var1    Var2    cor p   fdr
Lactobacillus_rogosae   Bacteroides_galacturonicus  0.467374392060265   0   0
Lachnospira_pectinoschiza   Bacteroides_galacturonicus  0.426779994497013   0   0
Lachnospira_pectinoschiza   Lactobacillus_rogosae   0.69419393866076    0   0
Bacteroides_galacturonicus  Bacteroides_salyersiae  0.074581964109831   0.548   0.829267117726658
Lactobacillus_rogosae   Bacteroides_salyersiae  0.320636458243079   0.002   0.03622024291498
Lachnospira_pectinoschiza   Bacteroides_salyersiae  0.132857736776534   0.21    0.56910941475827
Bacteroides_galacturonicus  Holdemanella_biformis   0.19490130543352    0.094   0.394447279549719
Lactobacillus_rogosae   Holdemanella_biformis   0.14472768883205    0.174   0.527328455284553

```

These are two correlation matrices that I want to bind through `rbind()`. The problem is that the exact correlation is represented in two ways in those two dataframes.

A bacteria is found in `Var1` of **Dataframe:1** while its pair is found in `Var2` of **Dataframe:2** (**see first correlation**).

I want to make a column so that each pair of correlations can be represented in the same format. e.g.:

```
Var1:Var2   cor p   fdr
Lactobacillus_rogosae:Bacteroides_galacturonicus    0.967374392060265   0   0
Lactobacillus_rogosae:Bacteroides_galacturonicus    0.467374392060265   0   0

```

Can someone help me achieve this, please?

# Answer
I think this should be fairly straightforward:

```
# working on Dataframe2
Dataframe2$tmp <- Dataframe2$Var2   # Create a temporary column of Var2
Dataframe2$Var2 <- Dataframe2$Var1  # Replace Var2 values by Var1
Dataframe2$Var1 <- Dataframe2$tmp   # Re-attribute tmp values to Var1
Dataframe2$tmp <- NULL              # delete tmp column

# Binding the two dataframes
df <- rbind(Dataframe1, Dataframe2)

# Concatenating the two columns
df$Var1 <- paste(df$Var1, df$Var2, sep = ":")
df$Var2 <- NULL                     # delete Var2 column

```

PS: I haven't tested this code so let me know if it works for you.

EDIT1:

```
InvertedRows <- Dataframe1[!(Dataframe1$Var1 %in% Dataframe2$Var2),]
InvertedRows$tmp <- InvertedRows$Var2   # Create a temporary column of Var2
InvertedRows$Var2 <- InvertedRows$Var1  # Replace Var2 values by Var1
InvertedRows$Var1 <- InvertedRows$tmp   # Re-attribute tmp values to Var1
InvertedRows$tmp <- NULL              # delete tmp column
Dataframe2<- rbind(InvertedRows, Dataframe1[Dataframe1$Var1 %in% Dataframe2$Var2,])
df <- rbind(Dataframe1, Dataframe2)

```
> 1 votes

---
Tags: r, correlation
---


# Question
Title: How do I find unique file names without their extensions in shell script?
I need to generate a list of unique file names without their extensions or using the shell script command "find"

Say I have a list such as:

```
file2022_col.R1.twx.gz
file2022_col.R2.twx.gz
file2022_col.bwa.bam
file2332_col.R1.twx.gz
file2332_col.bwa.bam
file3442_col.bwa.bam.bai

```

I want a list of unique file names such as:

```
file2022_col
file2332_col
file3442_col

```

I'm trying to use `find -type f -name` but I'm not sure how to construct the option to restrict the find command.

# Answer
It is simpler to use `cut` command to remove filename extensions:

```
ls -1 | cut -d '.' -f1 | sort | uniq

```
> 2 votes

# Answer
```
ls | perl -p -n -e  's/(.*)(\..*){1,3}$/$1/' | uniq

```

or

```
find . -name "file*" | perl -p -n -e  's/(.*)(\..*){1,3}$/$1/' | uniq

```

`find` *per se* is a really powerful command, but its a while since I've used it full functionality.
> 1 votes

# Answer
Use `find`, remove all extensions with `perl`, and make unique with `sort`:

```
find . -type f -printf '%P\n' | perl -pe 's{[.].*}{}' | sort -u

```

The Perl one-liner uses these command line flags:
`-e` : Tells Perl to look for code in-line, instead of in a file.
`-p` : Loop over the input one line at a time, assigning it to `$_` by default. Add `print $_` after each loop iteration.

`s{[.].*}{}` : replace `[.]` (a literal dot), followed by `.*` (any character repeated 0 or more times) with en empty string.

**SEE ALSO:**
`perldoc perlrun`: how to execute the Perl interpreter: command line switches
`perldoc perlre`: Perl regular expressions (regexes)
`perldoc perlre`: Perl regular expressions (regexes): Quantifiers; Character Classes and other Special Escapes; Assertions; Capture groups
`perldoc perlrequick`: Perl regular expressions quick start
> 1 votes

---
Tags: linux, shell
---


# Question
Title: How to make pan-core genome curve through command line on linux?
I´m working with a dataset of 566 genomes to analyze a pangenome. So I was working with PANWEB to create this pan core genome curve, however, there is too much sequence to work with this webserver. Well, specifically I´m looking for this kind of plot:

So, I just want to ask for some ideas of programs to do that with a big dataset and with a command line.

Here is the reference: https://www.researchgate.net/figure/Pan-genome-and-core-genome-Graph-representing-the-pan-genome-blue-and-core-genome\_fig1\_317390138

# Answer
> 3 votes
The command `create_pan_genome_plots` from Roary can create a plot like this:

---
Tags: genome, metagenome
---


# Question
Title: How does split reads look like in sam files?
I used bwa mem to align the DNA with the reference genome. If there are split reads (chimeric reads, come from two different parts of the genome), will it be split into two lines rather than one line?

# Answer
> 4 votes
I'm assuming you mean supplementary/chimeric alignment? The SAM Format Specification has a really detailed explanation as well as its Optional Fields Specification.

In general though, you will have an additional line in the SAM file to represent a chimeric alignment yes. For a single read, one line will represent the primary alignment; this is arbitrarily chosen. The other line will represent the split/supplementary alignment and have its supplementary alignment flag set. If there are multiple chimeric alignments, each one will have its own line.

You check if a primary alignment has any supplementary alignments associated to it by checking the SA tag which will list the location(s) of any supplementary alignment(s).

---
Tags: sequence-alignment, bwa, sam
---


# Question
Title: Can I Incorporate svaseq() into GSEA/GSVA analysis?
I understand GSEA/GSVA can take microarray-like expression matrix such as the output of voom() or vst(). However, I have a question on how we can also use svaseq() variables, correct the output from voom() or vst(), and apply GSVA on it.

I assume this is possible since I have read papers showing PCA after correction by SVA (such as this one and this one in their supplementary figure), but my statistical and coding knowledge does not allow me to recreate what they do. Is there a known pipeline that can do this

Edit: to clarify, I would like to generate a normalised matrix (e.g. rlog, or out put of voom from limma, or vst from DESeq2) that is adjusted for batch effect. ComBat, for example, can provide what I need. However, since I am doing my DESeq2 DGE analysis with svaseq(), I want to adjust it with svaseq variables. I am wondering if there are pipelines that can use svaseq vairables (or any varaibles) and adjust the value in the normalised expression matrix.

# Answer
Use `removeBatchEffects` from `limma`. The input counts should be on log scale, so `vst` and `rlog` are ok. The `covariate` argument can take the surrogate matrix from svaseq.

https://rdrr.io/bioc/limma/man/removeBatchEffect.html
> 1 votes

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


# Question
Title: Different results of spearman correlation between TPM and FPKM
TPM and FPKM of RNA-Seq data form GDC TCGA calculated based STAR were retrieved, respectively. The correlation between a specific gene, e.g. HIF1A, and other genes were calculated based on TPM and FPKM, respectively. And the significant genes were filtered based on P value \< 0.05. Here is the intersection of these two significant gene set shown in venn diagram. IMO, the two gene sets should not be obiviously different. Ang advice for this discrepance would be appreciated!

# Answer
This shouldn't be surprising that you see different correlations between gene expression data when expressed in different units. To see why, let's look at how these units are defined.

Let's denote the absolute number of fragments we sequence from an RNA-seq measurement as $X\_{gi}$ for gene $g$ in sample $i$, and $l\_g$ as the length of gene $g$ in kilobases. Then our different expression values will be:

$$ f\_{gi} = \frac{ X\_{gi} \cdot 10^6 }{ l\_g \cdot \sum\_h X\_{hi} } \text{ in FPKM} $$

or

$$ t\_{gi} = \frac{ X\_{gi} \cdot 10^6 }{ l\_g \cdot \sum\_h \frac{ X\_{hi} }{ l\_h } } \text{ in TPM} $$

You convert from FPKM to TPM units like so:

$$ t\_{gi} = 10^6 \cdot \frac{ f\_{gi} }{ \sum\_h f\_{hi} } $$

As you can see, *within* any single sample, $i$, this last denominator is fixed. But it is not fixed *across* samples (i.e. $\sum\_h f\_{hi} \ne \sum\_h f\_{hj}$ if $i \ne j$). This means that the ordering for a single gene's expression across samples is not preserved when you switch between units (i.e. if $f\_{gi} \> f\_{gj}$, it is **not** guaranteed that $t\_{gi} \> t\_{gj}$).

Because Spearman correlation is based on the order of each gene across samples in your dataset, and the order of genes in FPKM will be somewhat different from the order of genes in TPM, you will find different correlations between HIF1A and other genes when using the different units.

Here's a simplified, concrete example. Imagine that your transcriptome is made of up 3 genes (A, B, and C) with lengths 1 kb, 2 kb, and 3 kb, respectively. And imagine you're doing RNA-seq in 4 samples. If your absolute count matrix is this

$$ X = 10^6 \cdot \begin{bmatrix} 1 & 2 & 9 & 7 \\ 2 & 2 & 6 & 5 \\ 3 & 2 & 2 & 5 \end{bmatrix} $$

Then the expression matrix in FPKM is

$$ F = \begin{bmatrix} 166\\,666.7 & 333\\,333.3 & 529\\,411.76 & 411\\,764 \\ 166\\,666.7 & 166\\,666.7 & 176\\,470.59 & 147\\,058.82 \\ 166\\,666.7 & 111\\,111.1 & 39\\,215.69 & 98\\,039.22 \end{bmatrix} $$

and in TPM its

$$ T = \begin{bmatrix} 333\\,333.3 & 545\\,454.5 & 710\\,526.32 & 626\\,865.7 \\ 333\\,333.3 & 272\\,727.3 & 236\\,842.11 & 223\\,880.6 \\ 333\\,333.3 & 181\\,818.2 & 52\\,631.58 & 149\\,253.7 \end{bmatrix} $$

If you calculate the Spearman correlation for genes B and C against A, then you get these values.

| Gene | Correlation with A (FPKM) | Correlation with A (TPM) |
| --- | --- | --- |
| B | 0.316 | -0.8 |
| C | -1 | -1 |

You can see that gene C is anti-correlated with A, but B actually switches signs and goes from moderately correlated in FPKM to strongly anti-correlated in TPM.

So your observation that the same comparison on the different units gives pretty different results is expected. If you made a scatter plot with Spearman's $\rho$ from TPM along the x-axis and Spearman's $\rho$ from FPKM along the y-axis for all the genes that had a correlation p-value \< 0.05, you'll probably have a lot of genes that also flip signs like I showed with this toy example.

This paper by Wagner *et al.* is probably a good paper to review. It can give some context about RPKM and how to think about count-based sequencing data. Correlation-based analyses on data like this are typically used as a first-pass approach to get an idea about how things are behaving. But it's not very precise, as you've experienced. If you're doing exploratory work, I'd be inclined to pick TPM units, look at the results you have, then try to come up with a differential analysis question that can be handled in a more statistically robust way.
> 2 votes

---
Tags: rna-seq, correlation, fpkm, tpm
---


title: Grouping x-axis in scatter plot

text:
# Question
Title: Grouping x-axis in scatter plot
I am trying to plot one scatter plot for each group, where each group has continous values in x-axis. Reference Data

Gives continuous x-axis but all groups are combined.

```
ggplot(data_for_plotting,
aes(x = log10(Count), y = log10(Nonsilent.Mutation.Rate), color=factor(Clusters))) +
geom_point(alpha = 1, size = 1) + theme_bw()

```

Gives grouped x-axis but the variation in x-axis is lost.

```
ggplot(data_for_plotting,
aes(x = factor(Clusters), y = log10(Nonsilent.Mutation.Rate), color=factor(Clusters))) +
geom_point(alpha = 1, size = 1) + theme_bw()

```

Does achieve the grouping but the gives multiple plots.

```
ggplot(data_for_plotting,
aes(x = log10(Count),y = log10(Nonsilent.Mutation.Rate),color=factor(Clusters))) +
facet_wrap(~ factor(Clusters), shrink = TRUE) +
geom_point(alpha = 1, size = 1) + theme_bw()

```

I want x-axis to be grouped based on Cluster where each cluster has continous x-axis values.

# Answer
> 1 votes
You have two different types of values you want to display on the x-axis: `Cluster` (a categorical variable) and `Count` (a continuous variable). As Ram RS suggested, the easiest way to accomplish this is by taking your first plot (which uses a continuous x-axis scale) and facet along your discrete variable. This gets you what you want, and the rest of the work is adjusting the theme to make it look how you want.

I'm guessing based on your data that you're trying to create something that looks like this figure from Alexandrov *et al.*, Nature, 2013:

You can accomplish this by customizing the facet and the theme. Here's an example.

```
gg <- (
  ggplot(
    data = data_for_plotting,
    mapping = aes(
      x = log10(Count),
      y = log10(Nonsilent.Mutation.Rate),
      color = factor(Clusters)
    )
  )
  + geom_point(alpha = 1, size = 1)
  + facet_wrap(
    ~ Cluster,
    # this ensures that all facets are shown on the same row instead of wrapping
    nrow = 1
  )
  + theme_bw()
  + theme(
    # this removes the spacing between facets
    panel.spacing = unit(0, "lines"),
    # this removes the border around facets
    panel.border = element_blank()
  )
)

```

There are more parts of the theme you can adjust to your liking, like the facet header (e.g. `strip.text`). See ggplot2's reference documentation for details.

---
Tags: r, visualization, ggplot2
---


# Question
Title: Simulating 3' end tag-based scRNA-seq reads
Are there any tools that will simulate 3' end tag-based single-cell RNA-seq reads? For example, `Drop-seq`, `10X Chromium`, `CEL-seq`?

There are tools that simulate scRNA-seq gene count data (e.g. Splatter), but I can't find anything that will simulate reads or even counts at the transcript level.

# Answer
Looks like a couple options have become available:

**Minnow**: article code

**scReadSim**: preprint code
> 1 votes

---
Tags: scrnaseq, transcriptome, simulated-data
---


# Question
Title: How to convert 10x matrix.mtx.gz files to hdf5 format?
I have a bunch of folders containing barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz from the Cellranger Count output for a single-cell dataset that was sent to me from another lab. I need the .h5 version of these count matrices, but I do not have access to the original Cellranger Count output which contains the generated .h5 files.

Is there an R or Python function that can reliably convert these files into the proper hdf5 format? I haven't found anything online yet.

(x-post: https://www.reddit.com/r/bioinformatics/comments/veku5r/how\_to\_convert\_10x\_matrix\_tsv\_files\_to\_hdf5\_format/)

# Answer
Answered on biostars (https://www.biostars.org/p/9527636/#9527638) as well as on reddit: (https://www.reddit.com/r/bioinformatics/comments/veku5r/how\_to\_convert\_10x\_matrix\_tsv\_files\_to\_hdf5\_format/icqsfwn/)

The biostars answer is better as it needs just one package (DropletUtils) as opposed to the 2-3 required by the Seurat-allied packages suggested on reddit. Plus, Seurat can take forever to load a dataset.
> 0 votes

# Answer
Conversion into HDF5 from Python3 is straight forward, if this is a simple key: value dictionary format then it *should* be easy to mimic the CellRanger format. If it is a more complex data structure that would be tricky. The key to HDF5 is 'the HDF5 key' which you need to find (you would call it the directory)

To find the 'key' one approach is do is download https://www.hdfgroup.org/downloads/hdfview/ and then forward a screenshot of the CellRanger HDF5 file here. The directories are what HDF5 call 'keys' and you need these to write to and call from HDF5. There is a CLI for this operation.

An example of the data you are writing to would be needed.If per chance the data is in a Python dataframe you're in business, because these easily IO to HDF5

All you do is put the data into a simple dictionary and assign the dictionary to a key which is written to HDF5. Parsing the data from e.g. an AlignIO/SeqIO object to the dictionary is more complicated. When the HDF5 is read, the key has to be supplied then the dictionary is extracted and assembled back into an object.
> 0 votes

---
Tags: single-cell, cellranger, hdf5
---


# Question
Title: What is the contributing factor in single measurements of nanopore sequencing? One base, a k-mer or difference between the base left and entered?
My question is about nanopore sequencing and specifically about the current that is measured by the device in each measurement. The question is: In each measurement in nanopore sequencing, the change in the current between two sides of the membrane at the sensing region is recorded by the sequencer. Is it a single nucleotide or a k-mer or the difference between the first and last base in the kmer that contributes to one single measurement? I am citing two papers here that I looked at:

1. Nanopore sequencing technology, bioinformatics and applications (Nature 2021)
2. Three decades of nanopore sequencing (Nature 2016)

Based on my finding in these papers, I realized that all the nucleotide molecules in the k-mer in the sensing region contribute to the measurement, and not only one base, yet I am not sure if it is a good assumption.

I know some models like HMM and NN sequential models can computationally infer the individual bases (for example some basecallers). But my question is more about what is *physically* measured? In the second cited paper, it is stated that, in more recent nanopore designs, DNA strand can be fixed such that one individual nucleotide is positioned in the sensing region. Still, the diameter of the sensing region is not exactly one base. Moreover, since the current *change* is recorded not the current itself, the contributor could be the *difference* between the base which just passed from the sensing region and the one base that just entered?

So I restate my question; what is the contributing factor in single measurements of nanopore sequencing? In case it is the k-mer that is measured, how could we find k? Is it based on the properties of the device and the kit information?

# Answer
> 3 votes
The sequencing channels measure electric current, i.e. the flow of charged particles through the nanopore. The circuits are not measuring the bases directly; they're measuring the way the flow of ions is changed when different things are introduced into the pores.

Older R9 pores had a length that could fit about five bases in the pore at any one time, the newer R10 pores have a double-size pore that fits about ten bases. Any of these bases can potentially influence the flow of current, and that degree of influence is increased at the narrower pinch points of the Nanopore. Inconsistent unwinding of double-stranded DNA (or hybrid DNA/RNA) can also affect the current profile over time, so there may also be some effects on the current that transfer from where the unwinding occurs, 20 or so bases prior to the pore introduction.

The 5/6 kmer models are (or were, with R9) a reasonable sequencing model, because that's the number of bases that would fit in the pore, so those are the bases that have the greatest influence on the current.

I'll explain using a bigger physical model of hydroelectricity:

The mental model that I have in my head of a flow cell is a dam with a lake / reservoir, and holes in the dam that slowly drain the lake, spinning turbines that produce electricity (i.e. current). The interesting bits in the lake are the things going through the holes that aren’t water, but the water flow is the only thing that can be measured: there are no cameras attached to the holes, only current meters attached to the turbines.

It just so happens that those non-water things have a strange enough shape that they interrupt the water flow in large and predictable ways, which means that observing the water flow can be used as a good proxy for observing the things themselves.

If the holes get blocked, nothing can move through the holes, so there is no current flow: sequencing cannot happen.

Eventually that water reservoir will get depleted. If the reservoir is fully drained \[ionic balance of the carrier ions in the flow cell is equalised\], then it doesn’t matter how clean the holes are; the dam will produce no more electricity \[sequencing cannot happen\].

As an aside, this model helps me to understand why keeping the pores occupied with long DNA prolongs the sequencing life of the flow cell. The flow of carrier ions is reduced when a pore is actively sequencing, so the carrier ion reservoir drains slower, and can keep sequencing for longer.

---
Tags: nanopore
---


# Question
Title: the following errors are found during the installation of an open source software, UNCeqR on z shell
I am trying to do somatic mutation detection via the open source software, UNCeqR

I got following error was detected on mac z shell during installation.

```
BEGIN failed--compilation aborted at test.pl line 11.
not ok 1
make: *** [test_dynamic] Error 2
  CALLAHAN/Math-CDF-0.1.tar.gz
  /usr/bin/make test -- NOT OK
//hint// to see the cpan-testers results for installing this module, try:
  reports CALLAHAN/Math-CDF-0.1.tar.gz

```

I would like to ask for any advice to proceed.

# Answer
It is a Perl compilation error:

(Perl libraries) Math::CDF-0.1 -- Generate probabilities and quantiles from several statistical probability functions

> conda install -c bioconda perl-math-cdf

> cpan INSTALL Math::CDF

This module provides cumulative probabilities \[P(X \<= x)\] and quantiles \[given p, returns x such that P(X \<= x) = p\] from several statistical probability distributions.

This package provides a perl interface to DCDFLIB

https://metacpan.org/release/CALLAHAN/Math-CDF-0.1/source/README

`make` and `make test` of Perl libraries Math::CDF-0.1 in the installation are not okay

> https://drive.google.com/file/d/1Bd8TklsUPMqDRDn7zImGEXI9Nk84vkfU/view?usp=sharing
> 2 votes

---
Tags: software-installation, errors
---


title: beginner RNA-seq Replicate papers

text:
# Question
Title: beginner RNA-seq Replicate papers
have a good R and statistical analysis background (also with machine learning). in addition, I'm a fresh biotechnology grad. I would like to try to replicate some Rna-seq analysis with R papers (with their provided data). Any SHORT (beginner-friendly) papers to recommend?

# Answer
I like Gierlinksi, *et al.*, Bioinformatics, 2015 for this type of problem.

I like this paper because the methods are descriptive and clear, they have lots of figures that you can reproduce along the way, they spend a lot of time discussing "normalization" and quality control, and the data is easy to download.

Even better, the authors spent a lot of time making sure all the samples were split across lanes and flowcells in a consistent way so that both biological and technical differences could be robustly identified. They also compared a few different methods published at the time for performing differential expression analysis, so you can see how you can process the same data with different statistical models.

Unfortunately, they don't provide the explicit code they run for this paper, so you can't look at their code as a reference to guide you. But the methods are descriptive enough that I don't think there should be too much difficulty.
> 2 votes

# Answer
I am a beginner too. I will try the galaxy (training):

> https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/ref-based/tutorial.html
> 0 votes

# Answer
Here's https://drive.google.com/file/d/1Lwh6ofxi\_bb92mB0HoRrMJObE1MmtcMi/view?usp=sharing

It is from the book of `Methods of mathematical oncology`

https://github.com/ncbi/sra-tools

```

Suzuki, T., Poignard, C., Chaplain, M., & Quaranta, V. (2021). Methods of mathematical oncology : Fusion of mathematics and biology, osaka, japan, october 26-28, 2020 (1st ed. 2021 ed.). Singapore: Springer Nature Singapore : Imprint: Springer.

```
> 0 votes

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


# Question
Title: How to remove the unpaired reads in sam/bam files?
I have sam and bam files for the chimeric reads, which come from two different parts of the genome (For example, the first half of the read from part of Chromosome 1 and the second half of the read from part of Chromosome 3). I have removed the low mapping quality reads (MAPQ \<= 30) and reads with high mismatches. This might result in some unpaired reads which need to be removed before further analysis. May I ask whether there is any way to remove the unpaired reads in sam/bam files? Thanks!

# Answer
The `-f` option of `samtools view` is for flags and can be used to filter reads in bam/sam file matching certain criteria such as properly paired reads (`0x2`) :

```
samtools view -f 0x2 -b in.bam > out.bam

```
> 1 votes

---
Tags: sequence-alignment, samtools, phylogenetics
---


# Question
Title: How to align integer (non-DNA/protein) sequences?
I am looking for an algorithm to find the "best" alignment between two sequences of integers similar to how one aligns nucleic acids or amino acids for homology comparisons. For example, the best alignment for the two sequences below is:

```
                10 06 13 91 22 16 06
                |  |  |  |  |  |  |
       01 01 02 10 04 11 92 22 17 05

```

or an example with gaps:

```
                10 _  _  91 22 16 06 22 12
                |  |  |  |  |  |  |
          11 02 10 04 11 92 22 17 05

```

The alignment would be based on the distance between each number similar to aligning amino acid sequences where one would prefer to align a hydrophobic amino acid with another hydrophobic amino acid because they are more "similar". I hope this explanation makes sense.

Another way to think about this problem is if you took the 20 amino acids and replaced them with a number from 1-20 instead of their single letter abbreviation and then tried to align 2 sequenced based on how close the numbers where, how would someone align them? In this example the alignment is based on how close two integers are from one another whereas when aligning amino acids the alignment is based on the chemical similarity of the amino acids, ie leucine is very similar to isoleucine whereas valine is very dissimilar to tryptophan.

I've been experimenting with using the biopython pairwise library and trying to utilize it for my problem. There is a similar approach to what I am looking for here:

https://stackoverflow.com/a/45512865/1255817

This solution solves my problem, but only for the single character digits \[0-9\] and my numbers are more than one character.

Does anyone know of an approach (algorithm, package, library) to align sequences of integers? Ideally, I am looking for something in Python with a minimal working example.

# Answer
The tricky part about your request is the size of the alphabet. If you are happy with alphabets of up to about 250 "characters", then the programs in the FASTA package (github.com/wrpearson/fasta36) can be modified to do what you need. The are designed to work with unsigned char strings, and could work with matrices of up to 250 x 250 or so.

There is a non-vectorized version of Smith-Waterman that can calculate optimal local alignments with alphabets extended to 250 characters (with the appropriate scoring matrix).

However, I suspect that pairwise could do the same thing (with modifications), again using an expanded scoring matrix. The affine-gap Smith-Waterman algorithm is pretty simple, the trick is to get your multi-digit "characters" into a single byte (or unsigned char).
> 4 votes

# Answer
To parse data into `pairwise2` so you can follow the SO post here fairly easy by simply converting your output values into a higher base.

```
from numpy import base_repr
base36 = []
for count in range(1,20):
    base36.append( base_repr(count, 36))
myseq36i = str(''.join(base36))

```

*Output*

> 123456789ABCDEFGHIJ

This represents 1 to 20 as single digits (of course there's 36 single digit integers) and then implement the pairwise2 solution. The advantage of this solution is you're not stuck on 20 amino acids. `str` may not be needed.

The custom matrix described in that post you would supply can be based around the exact base 10 differences under the scoring matrix - to avoid subtracting in the higher base.

---

Please do upvote the answer on SO as it is relevant to your problem. The first solution they presented is cool and am sure there's broader application.
> 3 votes

---
Tags: phylogenetics, sequence-alignment, phylogeny
---


title: deseq2 makeContrasts function

text:
# Question
Title: deseq2 makeContrasts function
This article talked about various deseq2 designs etc. One of the designs I would like to use is explained as this:

Control versus treatment average

```
makeContrasts((treatmentI+treatmentII+treatmentIII)/3-treatmentCTL,
  levels=colnames(design))

```

I'm not sure how to incorporate this in deseq2 run which step?

Is it after running this

```
dds <- DESeq(dds)

```

Can I call the `makeContrasts` function?

# Answer
> 3 votes
The `makeContrasts` function is from edgeR so it does not help here. If you want to compare one group vs the average of several others here is a suggestion of how this could be done. It uses a design without intercept so all groups are easily accessable as coefficients in `resultsNames(dds)`:

```
suppressMessages(library(DESeq2))

#/ Example with four groups
set.seed(1)
dds <- makeExampleDESeqDataSet(m=12)
dds$condition <- factor(as.character(rep(LETTERS[1:4], each=3)))

#/ Use a design without intercept so all four groups are accessable as coefficients
design(dds) <- ~0+condition
dds <- DESeq(dds)
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing

#/ results using listValues to get A vs the average ov B,C,D
res <- results(dds, contrast=list(c("conditionA"), c("conditionB", "conditionC", "conditionD")), listValues=c(1, -1/4))

Created on 2022-06-21 by the reprex package (v2.0.1)

```

Determining whether or not it is appropriate to do this in a particular situation is up to you. You have to know what question you want to answer for your work. I cannot help with that.

---
Tags: deseq2
---


# Question
Title: BWA-mem and sambamba read group line error
*this question has been asked \[and answered\] on Stack Overflow*

This is a two-part question:

1. help interpreting an error;
2. help with coding.

---

1. I'm trying to run `bwa-mem` and `sambamba` to align raw reads to a reference genome and to sort by position. These are the commands I'm using:

```
bwa mem  \
  -K 100000000 -v 3 -t 6 -Y \
  -R '\@S200031047L1C001R002\S*[1-2]' \
  /path/to/reference/GCF_009858895.2_ASM985889v3_genomic.fna \
  /path/to/raw-fastq/'S\S[^_]*_L01_[0-9]+-[0-9]+'_1.fq.gz \
  /path/to/raw-fastq/'S\S[^_]*_L01_[0-9]+-[0-9]+'_2.fq.gz | \
/path/to/genomics/sambamba-0.8.2 view -S -f bam \
  /dev/stdin | \
/path/to/genomics/sambamba-0.8.2 sort  \
  /dev/stdin \
  --out host_removal/${SAMPLE}/${SAMPLE}.hybrid.sorted.bam

```

This is the error message I'm getting: `[E::bwa_set_rg] the read group line is not started with @RG`.

My sequences were generated with an MGI sequencer and the readgroups are identified like this `@S200031047L1C001R0020000243/1`, i.e., they don't begin with an @RG. How can I specify to sambamba that my readgroups start with @S and not @RG?

---

2. The commands written above are a published pipeline I'm modifying for my own research. However, among several changes, I'm not confident on how to define sample id as such stated in the last line of the code: `--out host_removal/${SAMPLE}/${SAMPLE}.hybrid.sorted.bam` (I'm referring to ${SAMPLE}). Any insights?

---

Thank you very much!

# Answer
Read groups must start with RG, that is a hardcoded requirement of the SAM format, nothing you can do about that.

See https://samtools.github.io/hts-specs/SAMv1.pdf
> 3 votes

# Answer
> My sequences were generated with an MGI sequencer and the readgroups are identified like this `@S200031047L1C001R0020000243/1`

Those are MGI read **names**, not read groups. (And the `@` is the FASTQ syntax that introduces a read name in a FASTQ file, not part of the read name or qname itself.)

Read **groups** more or less correspond to the sample that was sequenced. If a sample was sequenced across multiple sequencing machine runs and the resulting data merged into one file, there might be several read groups corresponding to the individual runs all pertaining to the same sample — though with increasing data volumes from individual sequencing runs this is perhaps less common than it used to be.

You get to make up the ids yourself, usually corresponding to what the sample is known as in your workflows. So your `-R` option should be something like

```
-R '@RG\tID:S123\tSM:S123\tDS:Sample 123\tPL:DNBSEQ'

```

(BWA's `-R` syntax is unfortunately not very easy to use. See addreplacerg for a similar `-r` option that is rather more user-friendly.)

The name you use for `${SAMPLE}` would usually be the same as the value of the `SM:` field within that `-R` option — the identifier that is used in your LIMS and workflows to identify that particular sample.
> 1 votes

---
Tags: sequence-alignment, ngs, genomics, bwa, sambamba
---


title: Normalization of data with RPkM

text:
# Question
Title: Normalization of data with RPkM
I'm having difficulty normalizing my data.

I was searching for transposable elements in my genome, and after this step, I made counts of reads in some transcripts. I produced something like this:

```
head(table_tissues_filtered_TE)
                  Lengths ova testes lobe retina suckers brain1 brain2 skin stage
Simple_repeat_80      134  58     77   48     69     115     137  131  195     75
tRNA_1                 59   0     14   12      1      19      12   14   21    104
Simple_repeat_87       26   1     33   12      3      15      24   21   19    180
Simple_repeat_114      22   0      0    0      1       0       0    0    2      7
Simple_repeat_115      30   0      0    0      0       0       0    0    0      1
Simple_repeat_123      22   2      3   317     45      13    652  651   15     21
                              axial                    gland viscera
Simple_repeat_80                 99                       35     557
tRNA_1                            9                        0       3
Simple_repeat_87                  9                        0       4
Simple_repeat_114                 0                      204       0
Simple_repeat_115                 0                       42       0
Simple_repeat_123               333                        5       4

```

where Lengths are the Length of each elements (simple repeats, etc), and the other columns indicate the reads counted with FeatureCounts. I have another table with the number of reads for each tissue:

```
head(reads_table)
       ova   testes lobe   retina  suckers   brain1  brain2  skin   stage   axial
      522444 310243 226146  102307  126055   489389  668243  372728 262536  233754
  gland  viscera
  24817   25689

```

I would like to do an RPkM analysis to normalize the data using R, but I don't know exactly how to do it. Can anyone can help me? thank you!!!

# Answer
> 1 votes
It's very unlikely that "a RPKM analysis" is the right answer. Assuming you'd like to do differential expression, using tools like DESeq or EdgeR on the count table are likely to be a better thing to do.

For reasons why RPKM is not a good approach, have a read of this answer.

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


title: ChrQ definition

text:
# Question
Title: ChrQ definition
What is ChrQ?

I assume it then means chromosme Q or something like that yet when I google for it I get very weird answers, but as I saw that a couple times when looking on freebayes pipelines manuals I wanted to be sure.

# Answer
*Answer from @devon-ryan and @terdon, converted from comments:*

Are you sure it's not just a made up chromosome name that's useful in examples since people won't blindly copy/paste the command and have it work?

It's almost certainly just a place holder, as Devon suggests, but we could give a more definitive answer if you could link us to the places you've seen this used. If you only see it in relation to freebayes then yes the author is using it as “foo”.

*\[please update / edit to improve this answer if you can\]*
> 1 votes

---
Tags: chromosome-capture
---


# Question
Title: How can extract the list of genes name from the raw data in GEO?
The format of the raw data that upload in GEO is different, like TXT, CEL or other form. How can I extract the **list of genes name** from these formats? And then how can I normalize them with R?

I would be really glad if anyone write **step by step** a complete set of codes for R with an example.

# Answer
> 0 votes
*Answer from @benn, converted from comments:*

You can find a step by step description in the limma user's guide. Just download the CEL files, and follow the user guide.

*\[please edit to improve this answer, if possible\]*

---
Tags: r, microarray, phylogenetics, data-download
---


# Question
Title: length of 'dimnames' [1686] must match that of 'dims' [3]
Please if anyone has experience with the use of the BSEQ-SC package for the deconvolution of bulk RNA sequencing data with single cell RNA sequencing data I will be very grateful for your suggestion.

I wish to deconvolve my bulk RNA seq data with single cell RNA seq data from the same tissue. I run in to the following error when running the code line below and I do understand the error but not how I can overcome it. Can any one please provide any useful hints?

```
csfit <- bseqsc_csdiff(ilc_bulk[genes, ] ~ treatment | fibroblast + macrophage,
                       verbose = 2, nperms = 100, .rng = 12345)

```

Then comes the error:

```
   Groups: a_nacl1=4L | a_nacl2=4L | a_nacl3=4L | a_nacl4=4L | a_nacl5=4L | NA0L
Cell type(s): 'fibroblast', 'macrophage' (2 total)
Fitting mode: auto
Data (filtered): 1684 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ...  OK
 Computing FDR using 100 permutations ... 101/100
  Alternative 'two.sided' ...   OK
  Alternative 'greater' ...   OK
  Alternative 'less' ...   OK
 OK
Timing:
   user  system elapsed
  2.136   0.153   2.267
Error in array(NA_real_, dim = c(nrow(b), ncol(cc), length(lev)), dimnames = c(rownames(b),  :
  length of 'dimnames' [1687] must match that of 'dims' [3]

```

```

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] edgeR_3.26.6                limma_3.40.6                forcats_0.4.0               stringr_1.4.0
 [5] purrr_0.3.2                 readr_1.3.1                 tidyr_0.8.3                 tibble_2.1.3
 [9] tidyverse_1.2.1             SoupX_0.3.1                 ggplot2_3.2.0               dplyr_0.8.3
[13] DropletUtils_1.4.3          SingleCellExperiment_1.6.0  SummarizedExperiment_1.14.0 DelayedArray_0.10.0
[17] BiocParallel_1.18.0         matrixStats_0.54.0          GenomicRanges_1.36.0        GenomeInfoDb_1.20.0
[21] RColorBrewer_1.1-2          xbioc_0.1.17                AnnotationDbi_1.46.0        IRanges_2.18.1
[25] S4Vectors_0.22.0            BisqueRNA_1.0               Seurat_3.0.2                preprocessCore_1.46.0
[29] e1071_1.7-2                 bseqsc_1.0                  csSAM_1.4                   Rcpp_1.0.2
[33] openxlsx_4.1.0.1            Biobase_2.44.0              BiocGenerics_0.30.0

loaded via a namespace (and not attached):
  [1] reticulate_1.13        R.utils_2.9.0          tidyselect_0.2.5       RSQLite_2.1.2          htmlwidgets_1.3
  [6] grid_3.6.1             Rtsne_0.15             devtools_2.1.0         munsell_0.5.0          codetools_0.2-16
 [11] ica_1.0-2              future_1.14.0          withr_2.1.2            colorspace_1.4-1       rstudioapi_0.10
 [16] ROCR_1.0-7             gbRd_0.4-11            listenv_0.7.0          NMF_0.22               Rdpack_0.11-0
 [21] labeling_0.3           GenomeInfoDbData_1.2.1 bit64_0.9-7            rhdf5_2.28.0           rprojroot_1.3-2
 [26] vctrs_0.2.0            generics_0.0.2         R6_2.4.0               doParallel_1.0.14      rsvd_1.0.2
 [31] locfit_1.5-9.1         bitops_1.0-6           assertthat_0.2.1       SDMTools_1.1-221.1     scales_1.0.0
 [36] gtable_0.3.0           npsurv_0.4-0           globals_0.12.4         processx_3.4.1         rlang_0.4.0
 [41] zeallot_0.1.0          splines_3.6.1          lazyeval_0.2.2         broom_0.5.2            modelr_0.1.4
 [46] BiocManager_1.30.4     yaml_2.2.0             reshape2_1.4.3         backports_1.1.4        tools_3.6.1
 [51] usethis_1.5.1          gridBase_0.4-7         gplots_3.0.1.1         sessioninfo_1.1.1      ggridges_0.5.1
 [56] plyr_1.8.4             zlibbioc_1.30.0        RCurl_1.95-4.12        ps_1.3.0               prettyunits_1.0.2
 [61] pbapply_1.4-1          viridis_0.5.1          cowplot_1.0.0          zoo_1.8-6              haven_2.1.1
 [66] ggrepel_0.8.1          cluster_2.1.0          fs_1.3.1               magrittr_1.5           data.table_1.12.2
 [71] lmtest_0.9-37          RANN_2.6.1             fitdistrplus_1.0-14    pkgload_1.0.2          hms_0.5.0
 [76] lsei_1.2-0             xtable_1.8-4           readxl_1.3.1           gridExtra_2.3          testthat_2.2.1
 [81] compiler_3.6.1         KernSmooth_2.23-15     crayon_1.3.4           R.oo_1.22.0            htmltools_0.3.6
 [86] Formula_1.2-3          lubridate_1.7.4        DBI_1.0.0              MASS_7.3-51.4          Matrix_1.2-17
 [91] cli_1.1.0              R.methodsS3_1.7.1      gdata_2.18.0           metap_1.1              igraph_1.2.4.1
 [96] pkgconfig_2.0.2        registry_0.5-1         plotly_4.9.0           xml2_1.2.1             foreach_1.4.7
[101] dqrng_0.2.1            rngtools_1.4           pkgmaker_0.28          XVector_0.24.0         rvest_0.3.4
[106] bibtex_0.4.2           callr_3.3.1            digest_0.6.20          sctransform_0.2.0      tsne_0.1-3
[111] cellranger_1.1.0       dendextend_1.12.0      curl_4.0               gtools_3.8.1           nlme_3.1-141
[116] jsonlite_1.6           Rhdf5lib_1.6.0         desc_1.2.0             viridisLite_0.3.0      pillar_1.4.2
[121] lattice_0.20-38        httr_1.4.0             pkgbuild_1.0.3         survival_2.44-1.1      glue_1.3.1
[126] remotes_2.1.0          zip_2.0.3              png_0.1-7              iterators_1.0.12       bit_1.1-14
[131] class_7.3-15           stringi_1.4.3          HDF5Array_1.12.1       blob_1.2.0             caTools_1.17.1.2
[136] memoise_1.1.0          irlba_2.3.3            future.apply_1.3.0     ape_5.3

```

# Answer
> 0 votes
*Answer from @stupidwolf, @benn, @llrs, converted from comments:*

The error complains about the 1686 rownames, while you have 1684 features. I think something is wrong with the naming of your rows.

This may be a bug that you should report to the package maintainer. I cannot really tell without more information about your eset. Looking at this line "Cell type(s): 'fibroblast' (1 total)", is it possible you only have one group under fibroblast? It will not make sense to run a random effect for only 1 group.

*\[please edit to improve this answer, if possible\]*

---
Tags: r, rna-seq, single-cell
---


# Question
Title: Elegant Way to Extract Gene Symbols from FASTA Headers
I've downloaded the RefSeq transcript sequences (FASTA format) for the GRC38 patch 14 build of the human genome and would now like to extract the gene symbols (only) in order from the headers. I'm confident I can do this using `grep` and regular expressions--the trick will be coming up with an appropriate pattern. (The gene symbol *appears* to be surrounded by the last pair of parentheses in a header, but I'm not sure this is absolutely true.) However, I'm wondering if there is a more elegant way to do this--e.g., via the NCBI EDirect utilities or `SeqKit`--that I'm overlooking. (I have a list of the accession.version numbers.) Any thoughts on this would be greatly appreciated.

# Answer
We can't be sure without seeing your file, but if what you really want is the text between the last parentheses on lines that start with `>`, then you can do this with GNU `grep` (the default on Linux) with:

```
grep -oP '^>.*\(\K[^)]+' file.fa

```

The `-o` means "print only the matching portion of the line" and the `-P` activates Perl Compatible Regular Expressions (PCREs). The regex looks for a `>` at the beginning of the line (`^>`), then as many characters as possible until the last parenthesis (`.*\(`), and then forgets everything matched up to this point with the `\K`. This means that the `-o` will only see things *after* the `\K`, a PCRE feature. So finally, we look for the longest stretch of one or more non`-)` characters.

If you don't have GNU `grep`, you can do the same basic thing with `perl` or `sed`:

```
sed -n 's/^>.*(\([^)]*\).*/\1/p' file.fa
perl -lne '/^>.*\(([^)]+)/ && print $1' file.fa

```
> 1 votes

# Answer
There are probably better ways of doing this with Python or another tool, but here's a quick method that works in bash (tested on Linux) assuming that the Gene Symbols are contained in parentheses, i.e. "(Gene Symbol)":

```
pat='\(([^)]+)\)'
while read line; do
    if [[ "$line" =~ $pat ]]; then
        echo ${BASH_REMATCH[1]};
    fi;
done < human.1.rna.fna

```
> 0 votes

---
Tags: fasta, phylogenetics, grep
---


# Question
Title: module differential connectivity (MDC) analysis
In this paper I am citing the lines below:

> Intra-modular connectivity of resilient and susceptible network modules was compared to that of corresponding genes in control mice to identify gain, loss, or no change.Pie charts (lower panel) summarize module differential connectivity (MDC) analysis. Proportionally, very few resilient modules had significant MDC compared to more than half of all susceptible modules, which predominantly showed gain of connectivity
>
> A module differential connectivity (MDC) metric quantified differences in coexpression networks to identify phenotype-specific modules (Zhang et al., 2013)

Need to clarify whether this "module differential connectivity (MDC)" analysis is a part of WGCNA or its an independent analysis. Because if I do a search for `module differential connectivity(MDC)analysis` the result I get is this \`DGCA: Differential Gene Correlation Analysis'. But no where in the paper they have mentioned "DGCA" library in the citation.

It would be helpful if i get some suggestion regarding this analysis.

This is the paper:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3677161/

# Answer
*Answer from @devon-ryan, converted from comments:*

They described how they did it in that paper, there's just not a tool for it. It's not part of anything. They describe the algorithm, maybe DGCA implemented it, maybe not, but they don't mention using it. WGCNA was obviously used to define the modules, though.
> 1 votes

---
Tags: r, wgcna
---


title: Construct phylogeny from a Fasta file

text:
# Question
Title: Construct phylogeny from a Fasta file
I have a set 189 taste receptor protein sequences (not aligned) in a fasta file. I like to get the phylogenetic tree in newick format. I was using earlier https://ngphylogeny.fr/. Unfortunately now it is not working due to some bugs. Necessarily, I just need newick information for the input fasta file.
Any help would be highly appreciated.

Here is one newick example:

```
((raccoon:19.19959,bear:6.80041):0.84600,
  ((sea_lion:11.99700, seal:12.00300):7.52973,
     ((monkey:100.85930,cat:47.14069):20.59201,
       weasel:18.87953):2.09460):3.87382,
  dog:25.46154);

```

# Answer
You will first need to align your sequences, and then feed the alignment to a program that can generate trees. I am very much out of date in this field and I haven't done anything like that in at least ten years, but back when I used to do it, my workflow was:

1. Align the sequences using MAFFT or t-coffee.
2. Open the alignment in the Jalview alignment editor/viewer which has options to generate the tree:

then just save the tree in Newick format:

The screenshots were taken using the default alignment provided by JalView itself, and this is the tree I got when I saved as Newick:

```
$ cat tree
((FER1_MESCR:24.484375,(((O80429_MAIZE:23.699999999999996,FER1_MAIZE:17.300000000000004):21.240625,(FER2_ARATH:13.067708333333336,((Q93Z60_ARATH:8.0,FER1_ARATH:8.0):8.401785714285715,(FER_BRANA:14.796875,FER3_RAPSA:8.203125):1.5982142857142847):0.9322916666666643):18.259375):9.75390625,((Q7XA98_TRIPR:24.8515625,FER1_PEA:27.1484375):6.643229166666671,((Q93XJ9_SOLTU:10.125,FER1_SOLLC:8.875):5.61363636363636,(FER_CAPAN:5.0,FER_CAPAA:0.0):19.38636363636364):20.48177083333333):8.97265625):9.52734375):20.2578125,FER1_SPIOL:20.2578125);

```
> 5 votes

# Answer
Here's the result from FigTree. I thought there was a bug in the newick file, but it was fine.

I strongly recommend the answer by @terdon because there are other issues here, for example `111.16` value for the root and the `100.86` which occurs twice isn't cool.

I should have placed `monkey` as the outgroup, anyway just please ensure you do this for the report.
> 1 votes

---
Tags: phylogenetics
---


# Question
Title: rRNA genes missing from metagenomic bins; is there a way to recruit them?
I have a couple of metagenomic bins which are okay in quality, but often missing rRNAs (16S, 23S...). I assume this is due to high population variability, combined with the high conservation of rRNAs, and reads of at most 300 bp. My binning software of choice is vamb.

Could this be fixed? Once a metagenomic bin is built, is there a way to see which 16S (or multiple 16S possible sequences) are likely to be associated with it? I considered using tools which explore the assembly graph, either manually (Bandage) or automatically (spacegraphcats), but I'm not sure if this is the right approach. After all, if the rRNA genes were linked in the assembly graph to genes that now clearly belong to a bin, they would have been added to that bin...It feels like this is a problem others may have encountered before, so I was wondering if anyone has an idea : )

# Answer
> 2 votes
I considered mentioning this but considered it implausible. However, apparently people are actually able to use PE shotgun data to link rRNA genes to MAGs.

This is something that is quite trivial with linked-read or Hi-C data, but I am assuming that you don't have these data (and don't have the desire to go back and generate them). However you presumably have PE shotgun data that you used to generate your MAGs.

---
Tags: rna, metagenome
---


# Question
Title: deseq2 model design : Different gene output
Here I have two designs:

**Design 1**

```
 dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered,
           colData=coldata, design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB)

```

To see the results I did this

```
resultsNames(dds)
 [1] "Intercept"          "Sex_Male_vs_Female" "Age"                "TMB"                "WBC"                "BM_percentage"      "FAB_M1_vs_M0"
 [8] "FAB_M2_vs_M0"       "FAB_M3_vs_M0"       "FAB_M4_vs_M0"       "FAB_M5_vs_M0"

```

I extracted only the Age from the first design

**Design 2**

I again used the design using `Age` variable :

```
dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered,
          colData=coldata, design= ~ Age)

```

As an example I take this gene this output from the first design

```
  gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
  <chr>              <dbl>          <dbl>  <dbl> <dbl>    <dbl>      <dbl>
1 ENSG00000196188     97.4         -0.742  0.113 -6.58 4.57e-11 0.00000125

```

This from the second design:

```
  gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
  <chr>              <dbl>          <dbl>  <dbl> <dbl>    <dbl>      <dbl>
  ENSG00000196188     97.4         -0.757  0.119 -6.34 2.28e-10   3.43e-07

```

Another example this gene from the first design; the result is not significant:

```
           gene  baseMean log2FoldChange  lfcSE    stat pvalue   padj
ENSG00000173826  74.54467     0.01283834  0.589  0.0218  0.983  0.993

```

The same gene from **Design 2** *is* significant:

```
           gene  baseMean log2FoldChange  lfcSE    stat   pvalue     padj
ENSG00000173826 74.54467       -2.004743  0.190   -10.6 4.76e-26 1.35e-21

```

Is this drastic outcome is due to the maths behind it or am I doing something wrong?

So would my interpretation be as such when i consider Age as my only variable I can consider `ENSG00000173826`significant ? But it is not significant when my main effect which i want to test is FAB which is my design 1.

I would like to know how to interpret and report the same.

Any suggestion or help would be really really appreciated

# Answer
It's not surprising to me that the adjusted fold change statistic changes in different ways for different genes when covariates are added or removed. That's sort of the point of using covariates, as long as it makes sense to include them (which is an experimental decision, not a mathematical or statistical decision).

One thing I'd recommend checking is to make sure the results you're looking at are what you *think* you are looking at. You say that you're interested in testing `FAB`, but that is not included in the second (Age-only) design, which means that it won't be possible to extract any `FAB` comparisons from the results. It looks like your first design might be looking at `Sex` (based on the resultsNames ordering) or one of the `FAB` comparisons (based on reverse ordering of variables in the model), whereas the second might be looking at `Age`. This should be explicitly stated when displaying the results, but you haven't included that information in your output.

If you want `FAB` comparisons in the results, that variable needs to be in the specified model, e.g. "design = ~ Age + FAB". The Design 2 model that you have in your question only includes `Age`.

A good way to do this without getting the wrong contrast is to fully specify the variables when fetching results, as explained in the DESeq2 vignette:

> In any case, the `contrast` argument of the function results takes a character vector of length three: the name of the variable, the name of the factor level for the numerator of the log2 ratio, and the name of the factor level for the denominator. The contrast argument can also take other forms, as described in the help page for results and below

```
res.age <- results(dds, name="Age")
res.FAB.M3.M0 <- results(dds, contrast=c("FAB","M3","M0"))

```
> 2 votes

---
Tags: deseq2
---


# Question
Title: Filtering paired-end reads with sambamba: avoid discarding reads on the minus strand
I have a BAM file (DNA, shallow whole genome sequencing at ~1X) where I want to filter reads (using sambamba) to keep only those which have a template length \> 20 and mapping quality \> 20, discarding unmapped reads.

The filter I have used is this:

```
sambamba view -h \
   -t 5 \
   -F "paired and template_length > 0 and mapping_quality > 20 and not (unmapped or mate_is_unmapped)" \
   -f bam 142-ready.bam \
   -o 142-filter.bam

```

However, inspection of the reads in IGV led me to think that perhaps reads on the minus strand are being somehow filtered out.

Pre filter:

Post filter:

The only way I don't get such an aggressive filtering is by setting only `not (unmapped or mate_is_unmapped)`. But that, obviously, defeats the purpose as it doesn't include what I need.

I can't, however, see what exactly is wrong with the filter. What is the correct way to filter for template length while not accidentally discarding reads on the minus strand?

# Answer
> 1 votes
The reason is that `sambamba` does not use the absolute value of the template length, so, if the read maps so the minus strand, the template length is actually negative. Setting a template length \> 0 means that all reads mapping to the minus strand are removed.

In fact, you can't still filter for negative template lengths. See this GitHub PR.

---
Tags: sam, filtering, reads, sambamba
---


title: What is the function of the heteroatoms

text:
# Question
Title: What is the function of the heteroatoms
I am trying to make a prediction application using the sequence of proteins.

If we think about the 6COU protein, I'm running the code below

```
pdb = PDBParser().get_structure("6c0u", "6c0u.pdb")

for residue in pdb.get_residues():
    print(residue)

```

Some of the output looks like below when looking at residues. What is the function of the heteroatoms here, and follow up question;

When reading a paper on Binding Site Detection of Proteins, there is exactly one sentence in the paper.

***"We first clean the input structure by removing all heteroatoms and solvent molecules from the protein structure using the Biopython47 library."***

What is the meaning of this sentence if you are working on protein sequences, heteroatoms are not already there, what does it mean to delete them?

If you want to check the paper, here

# Answer
The text you cited mentions "input structure", so they work with structures, not just "protein sequences".

Heteroatoms = atoms marked as `HETATM` in the PDB and mmCIF formats. In these formats each atom is listed as either `ATOM` or `HETATM`. These are not particularly meaningful categories, but they are there, so they are easy to use. It's a common practice to remove heteroatoms with `grep -v ^HETATM`.

Note that in the files from the Protein Data Bank not only solvent and ligands are marked as HETATM, but also "non-standard" amino-acids in a protein. An amino acid is considered standard if it can be genetically encoded. So in the files from the PDB selenocysteine is ATOM and selenomethionine is HETATM. Some programs don't follow these rules and output all protein residues as ATOM.

Removing heteroatoms works well for the vast majority of structures that have no modified residues. But in some cases it will remove also some protein residues, which may not intended.

**Update**: if you are wondering **why** there are ATOM and HETATM records, I have never came across the rationale. The PDB format specification has had both types of records since 1976. Different versions of the format spec use different wording, but the meaning is always similar:

* ATOM records give atomic coordinates for "standard" groups.
* HETATM records give atomic coordinates for "non-standard" groups.

(Although standard residues that not in a polymer are also HETATM.)

In the mmCIF format the same record is preserved in \_atom\_site.group\_PDB with the description:

> This data item is provided for compatibility with the original Protein Data Bank format, and only for that purpose.

Most of the structures have only standard residues in protein/DNA/RNA polymers, so removing HETATM records removes only ligands and solvent.
> 6 votes

---
Tags: proteins, python, protein-structure, pdb
---


title: HDF5 and BioSQL solutions

text:
# Question
Title: HDF5 and BioSQL solutions
I'm looking at better database/storage solutions for NCBI virus data, with all attributes particularly year and country of isolation, together with structural data, possible antibody data, T-cell data and bioinformatics data such as %GC etc... The amount of data is comparatively small, with a maximum size of `3 x ([6000 < no. viruses < 20000] for 10 genes)` and each entry will contain nucleotide, amino acid and one or more pieces of alignment data (i.e. with indels).

The data will be used for numerous different calculations from structural biology, phylogenetics, and deep learning. Some of the calculations will require parallel processing. A typical calculation would be to identify a clade/similarity/difference then screen the database to look at the isolates information.

Previously I've used lots of flat files which are loaded at initiation and are held within a data pipeline either as an object (Python) or especially a multi-dimensional hash (Perl, preferred).

The two solutions that I think look suitable are:

1. BioSQL \- The official Genbank solution, however my SQL commands are very basic and the documentation appears minimal (for BioSQL)
2. HDF5 (hierarchical data format) - this looks good because it uses complex hashes (dictionaries), and allows layering of multiple information onto a single Genbank sequence code (key). It has strong documentation. However, I understand this isn't a "database" (long term storage), but thats how I intend to use it.

**JSON/Pickling** I had a look at JSON, which to me looks okay, because I'm fond of complex hashes / dictionaries, but it appears to be more of a Web solution. Finally, "pickling" (Python) is cool, very cool, but I understand this not a database solution, because it risk compatibility issues, e.g. with future Python releases.

I am interested in knowing what archival systems board members are using, and a critique whether my assessment is ok for moving forward.

# Answer
> 1 votes
Just to conclude on the solution to this ancestral question. HDF5 is very useful, but there are better solutions for this specific question if it is focused on pandas dataframes. If you've lots of data the pythonic way to do this is pandas and therefore storing as a pandas dataframe makes sense.

For long term storage, which was the original question, the storage of choice for pandas is `parquet`. Described here

> Apache Parquet provides a partitioned binary columnar serialization for data frames. It is designed to make reading and writing data frames efficient, and to make sharing data across data analysis languages easy. Parquet can use a variety of compression techniques to shrink the file size as much as possible while still maintaining good read performance.

In the following example, two methods `parquet` methods are used for a pandas dataframe proteinDF which is first written and then imported.

```
from fastparquet import write
write('outfile.parq', proteinDF)

```

and

```
import pandas as pd
pd.read_parquet('proteindata.parquet', engine='fastparquet')

```

The above approach appears to have received most current developmental attention. However,

```
import pyarrow as pa
import pyarrow.parquet as pq
arrow_table = pa.Table.from_pandas(proteinDF)
pq.write_table(arrow_table, 'myoutpath', use_dictionary=False,
                   compression='SNAPPY')

```

and

```
import pandas as pd
pd.read_parquet('example_pa.parquet', engine='pyarrow')

```

The above approach is compatible with `numba`. 'Multithreading' is further described here.

The differences between the methods appear slight. Parallelisation is a big deal and depends on the problem.

A general discussion on the methods is here

**HDF5 vs parquet** The backdrop for `parquet` is the method is used for batch I/O not record by record writing to disk and that is this question and for that is pretty fast (at a guess bit better than `pickle`, but a lot more compression). `HDF5` is perfect for rapid record by record writing and that is what it is designed for, so 'benchmarks' comparing I/O speed are not realistic unless they compare batch vs. record by record timings, where `HDF5` excels. `parquet` does achieve significantly better data compression and is considered a long term storage solution. In my hands `HDF5` for fasta writing results in around half the size of the flat file (not much).

# Answer
> 1 votes
Just to conclude, @user172818 and @AlexReynolds are collectively saying either learn SQL or stick with CSV.

---

Update, I can provide an answer. It's HDF5 with pandas dataframes (Python).

Simply

```
df.to_hdf('path/file.h5', 'key', format='table')

```

This isn't SQL but it works:

What it gives,

* it's fast, that's it's purpose (not my main reason though)
* it perfectly preserves all pandas dataframe information, this is important
* The 'key' ( object) allows you to store multiple dataframes into the same file, you just provide different key names

The last point is very useful, because you can pull out a dataframe of choice and the SQL commands are mimicked within pandas

I ran into trouble using CSV to store/recover complex pandas dataframe formats.

So it's good enough and much faster than SQL. I could have problems with format changes in time (unlike SQL), but for the here and now it works well.

---

Just to add a caveat, for dictionaries the json dump works well, but it is just shunting dictionaries around, so isn't relevant.

---
Tags: python, perl
---


# Question
Title: How to normalise the histogram height in matplotlib
I am trying to create a histogram plot with `matplotlib`. I am using `bmi` column from four different data frames (d1, d2, d2, d4).

The issue is that these 4 datasets have different sample sizes. So there is a disparity between the heights of the bar.

Is there any additional option in the matplotlib to normalize the heights of the bar?

import matplotlib.pyplot as plt

def bmiplot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/110.0, graphHeight/150.0), dpi=100)
    axes = f.add_subplot(111)

    d1 = dm_group1
    d2 = dm_group2
    d3 = dm_group3
    d4 = dm_group4

    plt.hist([d1.bmi.values.flatten(), d2.bmi.values.flatten(), d3.bmi.values.flatten(), d4.bmi.values.flatten()],
    label=['bmi_dm1', 'bmi_dm2', 'bmi_dm3', 'bmi_dm4'])
    plt.legend()
    plt.rc('xtick')
    plt.rc('ytick')

graphWidth = 700
graphHeight = 800
bmiplot(graphWidth, graphHeight)

```
</code></pre>

# Answer
One approach is `density`, traditionally the code would be ...

```
kwargs = dict(histtype='stepfilled', alpha=0.4, normed=True, bins=50)

plt.hist(d1, **kwargs)
plt.hist(d2, **kwargs)
plt.hist(d3, **kwargs);

```

`alpha` is transparency allowing overlapping distributions.

**Note** `normed` is deprecated, the blurb here describes `density` now replacing it.

The blurb is here:

> matplotlib.pyplot.hist(x, bins=None, range=None, density=False, weights=None, ... \**kwargs)

The definition of `density` is here

> densitybool, default: False If True, draw and return a probability density: each bin will display the bin's raw count divided by the total number of counts and the bin width (density = counts / (sum(counts) * np.diff(bins))), so that the area under the histogram integrates to 1 (np.sum(density * np.diff(bins)) == 1). If stacked is also True, the sum of the histograms is normalized to 1.

A different approach is relative frequency, via the `weights` command.

```
 plt.hist(mydata, weights= *)

```

\*, relative counts per data type / total counts over all data types
> 1 votes

---
Tags: python, graphs
---