# Question
Title: Approaches for filtering bacterial/fungal contaminant sequences from RADseq results
I'm working with RADseq data (288 compressed FASTQ files, `*.fq.gz`) from 3 plates of plant tissue. I've done demultiplexing and preliminary QC (using the Stacks `process_radtags` command), and before performing a reference alignment, I want make sure the sequences in each sample (roughly) correspond to the organism I expect them to be derived from, rather than potential contaminants (bacteria and/or fungi).

Does anyone have a straightforward, "canned" approach for doing this, either in R or a different environment? I believe the `blastn` algorithm can be used to query the sequences in these files against a database (like NCBI), and is usable on the command line. However, I've not used it before, and I'm not sure if it makes sense to query *all* the sequences in each file. Additionally, I don't know if the command can process this many files of this size (~600 MB per FASTQ file).

This post from ~4 years ago asks a similar question for a pair of human blood samples, but is concerned with bacterial/viral DNA. The answer linking `kmerfinder` might have promise, but I'd appreciate any input from folks familiar with this issue.

(*Background*: I've gotten my Master's degree in Biology and have worked with NGS data for ~4 years now. I'm performing analyses on a 32-core Linux server. I'm less skilled in Python, but am willing to use it if there are no other options.)

# Answer
I ultimately addressed this issue using the kraken2 software. I found this approach to rely on a good deal of my server's memory, but it was sufficient for doing a viral/bacterial classification of sequences.
> 1 votes

---
Tags: ngs, filtering, bacteria, blastn, radseq
---


# Question
Title: Converting convert .svs files into .mrxs format
I wish to convert Aperio AFI microscope .svs files from a pathology lab microscope into 3DHistech .mrxs format. Lab ppl have some pipeline X and they asked me to look into the matter.

I can read the SVS files via openslide or openmicroscopy. Conversion to 'regular' pictures (tiff, for instance). I figure out that since openslide can read mrxs, writing this format should be feasible too, but can't find any (free) software to do this. I don't care about any specialties, just need the end result to be readable by Histech readers.

Can anyone help with the conversion? Lacking full software, I could try to code writer if resulting file description is available.

# Answer
If you truly "just need the end result to be readable by Histech readers" the simplest (and maybe even cheapest in the long run) solution might be a paid ($) license for CaseViewer (or whatever their newer viewer is called). I have no idea what that would cost but many of these vendor wsi viewers can read and even write other formats for a price.
> 2 votes

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


# Question
Title: Are there any simple ways to call copy number variants from a Nebula Genomics WGS data?
I want to find my copy number variants. I got DNA testing from the Nebula Genomics and received the raw datas in .cram, .fastq, .vcf formats. Are there any simple and easy method to call my copy number variants from one of these files? I tried GATK, but I was unable to collect Read Counts from 1000 genomes datas.

# Answer
Yes and no. Yes, there are tools that can do this, but no, they do not work very well. Short read NGS sequencing just isn't the right tool to find large structural variants like CNVs. As a rule of thumb, you can expect around 30% accuracy when compared to other experimental approaches such as MLPA or array CHG. That said, yes, there are indeed various tools that can do this. My personal favorite is delly because it is accurate, actively maintained and since it is written in C it is quite fast.

Once you have downloaded and installed it, you'll likely need to convert your cram to bam as I don't think it can deal with cram:

```
samtools view -b foo.cram > foo.bam

```

Then you can run delly on it with (assuming your cram was aligned against hg38):

```
delly call -x path/to/delly/excludeTemplates/human.hg38.excl.tsv \
   -g path/to/hg38.fa -o out.bcf foo.bam

```

You will also need to download a copy of the human genome fasta file that Nebula used to generate your cram. You can ask them where to get one, and change `path/to/hg38.fa` with the path to that file.
> 0 votes

---
Tags: wgs, cnv
---


title: Population Allele Frequency

text:
# Question
Title: Population Allele Frequency
In UK biobank research environment

In GENOMICS tab of UK BRAP, if I type `7:140476600-140480852` coordinate in search bar and the same time after filtering my cohort for this coordinate, in `ALLELE TABLE` I see columns like:

```
Cohort AF Cohort Allele Count Population AF Population Allele Count GnomAD AF

```

I Googled a lot but no clue where population allele count and frequency come from

Does anyone have any experience with this?

# Answer
I don't know specifically about the UKBB RAP, however, Population Allele Count is the count of that specific allele genotyped in the population in question. In this case, it most likely refers to the population that your cohort is part of. The population allele frequency is the same thing but noted in terms of frequency, i.e. divided by 2 times the number of genotyped individuals in that population.
> 3 votes

---
Tags: allele-frequency, population-genetics, brap, dnanexus, biobank
---


# Question
Title: Does the `mmcif/pdb` format store atomic coordinates in units of Ångstrom, nanometers or something else?
I was surprised as to how hard it was to find a reference on this among the docs.

Does the Protein Data Bank mmcif/pdb-x archival format use Angstrom to specify cartesian x,y,z coordinates of atoms or something else or abstract units (perhaps defined per each `mmcif` `_data` frame?

How is the origin of the grid calculated -- for crystallographers out there?

Thank you

# Answer
> 4 votes
From google I got this as first hit :

PDB\_format\_1992.pdf :

Hopefully they keept the same Unit of measurement during the data/file format evolution/extension

It seems the case for `mmcif` Data Item \_atom\_site.Cartn\_x :

Some history here :

The Protein Data Bank: a historical perspective , The Protein Data Bank

AND :

```
PY  - 1971
DA  - 1971/10/01
TI  - Crystallography: Protein Data Bank
JO  - Nature New Biology
SP  - 223
EP  - 223
VL  - 233
IS  - 42
SN  - 2058-1092
UR  - https://doi.org/10.1038/newbio233223b0
DO  - 10.1038/newbio233223b0
ID  - ref1
ER  -

```

About the axis origin CHARMM-GUI has a nice introductory lesson to PDB format (see ***CRYTST1 Record*** Section) Lesson 1: Introduction to PDB Format

and stolen from the web, free part of an old book :

---
Tags: pdb, 3d-structure, distance
---


# Question
Title: Concatenate multiple fasta files into individual files per unique entry
I have several splitted fasta files like this:

```
TK273759.fastaaa
TK273759.fastaab
TK273759.fastaac
TK273759.fastaad
TK273759.fastaae

```

I want to concatenate these fasta files into one file like this:

```
TK273759.fasta

```

What command should I use to make this happen? I appreciate your help.

# Answer
> 0 votes
Please try searching online before asking for help. A simple `cat` will do the trick here:

```
cat TK273759.fasta[a-z][a-z] > TK273759.fasta

```

---
Tags: genome, bash
---


title: Defining out-groups in MEGA

text:
# Question
Title: Defining out-groups in MEGA
I am using the MEGA software for trimming protein sequences and building phylogenetic trees.

How can I define out-groups of phylogenetic tree in MEGA? Is there a way to calculate the distances of each species from the root?

# Answer
> 0 votes
Hmmm.... the issue is that the precise location of the root - unless it's a molecular clock - is arbitrarily assigned within limits. You can designate populations ('clades') and ask for the genetic distance between the clades, it's one of the menus header options. To do what you want you need to build the tree (neighbor-joining method) and assess it in the pop-up window. Switch the view of the tree to 'genetic distance' and then manually add the branch lengths for each root - tip distance. One of the Beast supporting packages use to calculate the root tip distance, I think this option was removed in the upgrades.

---
Tags: phylogenetics
---


title: LD local clump-GRCh38 bfile

text:
# Question
Title: LD local clump-GRCh38 bfile
I want to do LD local clumping. The tutorial (https://mrcieu.github.io/ieugwasr/articles/local\_ld.html) teaches how to, but the downloaded bfile is GRch37. I followed this question to download GRch38. LD clump GRCh38 GWAS results

I downloaded and it is "all\_hg38.pgen". How can I transform it into bfile? Thank you.

# Answer
I have the same problem when facing the analysis of eQTL data for clumping and I do not find any answer to solve it. So I try to find a method to transform 1000 Genomes phase 3 ailgned by hg38 data which is pfile format (pgen, pvar, psam). I found a method at: https://cran.r-project.org/web/packages/snpsettest/vignettes/reference\_1000Genomes.html. This is a method about how to extract the EUR population and transfrom to bfile format using PLINK 2.0. After this transfrom we can use PLINK 1.9 to clump the GWAS data. I think this is the method which we want to use and I will try it. Hope this can help you.
> 0 votes

---
Tags: r, ld
---


# Question
Title: Aligning sequences with multiple genetic codes!
I am doing a project on duplicated genes and I have a major difficulty on how to align sequences that use different genetic codes. I work with fasta files that contain sequences of protein coding genes, every fasta file includes genes from many species that are orthologs/paralogs. These genes even in the same species may use different genetic codes. My goal is to align the sequences, codon alignment \[the process of codon alignment considers the reading frame of the translated protein and adjusts the nucleotide alignment so that the protein alignment stays in frame.\] , and afterwards to create N-J trees. For the whole process I work on MEGA11. I cannot organize the species-genes into groups with common genetic codes because I need all of them in the tree. One option, to my mind at least, would be the alignment not in codon but in DNA,or another option is to align their protein sequences or lastly a third option could be the change of the sequences in order to get rid of the stop codons and afterwards to use one common code for everything, by doing that I know that I will lose some info but at least I can work.Do you have any idea how to overcome this obstacle? Do you agree with any of this options? Thank you in advance!

# Answer
> 0 votes
I am having a very similar issue, however mine are multiple sequences of the same specie. I have a mixture of DNA and RNA sequences so although it appears to aligning fine the proteins produced are not correct. I’m at a loss on how to fix this so if you managed to figure out what to do please let me know.

---
Tags: phylogenetics, sequence-alignment, multiple-sequence-alignment, codon, megax
---


# Question
Title: ipdSummary command stops working if --identify m5C_TET is included in command, and does not give any m5C info either
I have installed the latest SMRT Tools (ipdSummary version 3.0). My goal is to extract ipdRatio and modified bases. In the ipdSummary --help command :

```
--identify IDENTIFY   Specific modifications to identify (comma-separated list). Current options are m6A, m4C, m5C_TET. Using --control overrides this option. (default: m6A,m4C)

```

Accordingly, I am running this command:

```
ipdSummary sorted_chr21_000.bam --reference $REF_path/$REF_FILE_NAME.fasta --referenceWindow 20:1-200000 --identify m6A,m4C,m5C_TET -j 20 --methylFraction --outfile $save_path/$chrom/$name

```

My bam file is sorted and relatively small (\<200 MB). But this command never finishes and does not throw any errors. In different stages of working it just stops. I checked different thread numbers and smaller files, but still, the command doesn't work. However, once I changed 'm5C\_TET' to 'm5C' it worked perfectly fine and very fast on the same bam file, but in the output, I could not see any information on 'm5C' modifications. I checked the new command with many other bam files and still the same.

My question is the following: How can I get at least 'm5C' modifications? Should I change my SMRT Tools version?

# Answer
I asked the same question on the ipdSummary GitHub page and I got the following response:

"I think that was an experimental option that was never fully working or validated, and it no longer works. We probably need to update the help text. Use the separate tool jasmine if you want to do m5C calling."

jasmine also is a tool in the SMRT toolkit that analyzes the kinetic signatures of cytosine bases in CpG motifs to identify the presence of 5mC.
> 0 votes

---
Tags: pacbio
---


# Question
Title: Featurecount output .txt file from bam file
After running this command:

```
featureCounts -p -s 2 -a $genome -o $dir"/"$specie"/count_table/"$value".txt" $output_loc -T 8

```

on the BAM output of HISAT (`output_loc`), I got this output:

Now, does the last column of this file show the gene count value? And if I merge the same column of different samples, how should I go through the normalization step?

Thanks a lot for your help!

# Answer
The first 6 column in standard featureCounts output represent what is in the column names. All columns after than (starting at 7) represent the counts for the sample(s). If you use a single bam file as input then it's one column, if you use many bams as input then it is one column per bam. If the counts are gene level, exon transcript or anything else depends on how you run the tool. Typically, with default settings, it is usually the raw gene level counts.

Generally, you can sum up technical replicates (=lane/sequencing replicates to increase depth) prior to running normalization.
> 3 votes

---
Tags: rna-seq, gene-expression, read-mapping, featurecounts
---


# Question
Title: Printing the results of for loop into a txt file
I have a for loop which prints available vcf files in the path for a list of patient IDs **(using find dx data function)** into my zsh terminal:

```
for i in $(cat patient_id.txt); do
 dx find data --property patient_id=$i --path "vcf/"
done

```

How I can have a file contain the vcf files that the loop shows in terminal? I mean by using this command in the terminal, I am seeing the name and path of my vcf files for a given participant list and working directory; Instead of seeing the file name in the terminal, I want to have that file name in a txt file.

Thanks

# Answer
This command worked just by adding `> file` after `done`

```
for i in $(cat patient_id.txt); do
dx find data --property eid=$i --path 'vcf/'
done > vcf_list_file.txt

```
> 3 votes

# Answer
Things that write to standard output on the terminal can be wrapped up in brackets and redirected to a file. Assuming this is what you want to do, I added brackets around the existing code, then used `>` to redirect the output to a file.

```
(for i in $(cat patient_id.txt); do
 dx find data --property patient_id=$i --path 'vcf/'
done) > vcf_files.txt

```

Note that the brackets go around the entire code, including the `for` bit.
> -1 votes

---
Tags: bash, dnanexus, zsh
---


# Question
Title: Subsetting a Seurat object based on colnames
I'm trying to subset my seurat object based on colnames. I have gone ahead and labeled each cluster and now I want to subset all the colnames that are in Cancer\_human for human\_colnames and all the barcodes that are not in Cancer\_human for mouse\_colnames but get an error.

```
human_colnames = colnames(scData[,scData$cell_ann == "Cancer_human"])
mouse_colnames = colnames(scData[,!scData$cell_ann == "Cancer_human"])

Error in `[.Seurat`(scData, , scData$cell_ann == "Cancer_human") : Incorrect number of logical values provided to subset cells

```

I'm filtering on colnames not rownames. I want all the barcodes from the Cancer\_human.

# Answer
Did you tried the `subset` function?

```
# assuming that scData@meta.data looks like this
#         orig.ident  ...  cell_ann
# ATTCAA  sample1          Cancer_human
# ATTTAT  sample1          Cancer_human
# ...
# AAATCA  sampleX          Cancer_mouse
# ATTTAA  sampleX          Cancer_mouse

Idents(scData) <- "cell_ann"
human_subset <- subset(scData, idents = "Cancer_human")
mouse_subset <- subset(scData, idents = "Cancer_mouse")
mouse_subset <- subset(scData, idents = "Cancer_human", invert = TRUE) # alternative way

```

P.S. It is not clear what you are trying to do. You should add to your post how your data looks like.
> 0 votes

# Answer
If you want the column names to be preserved, you can subset the column names as a vector:

```
human_colnames = colnames(scData)[scData$cell_ann == "Cancer_human"]
mouse_colnames = colnames(scData)[!scData$cell_ann == "Cancer_human"]

```
> 0 votes

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


# Question
Title: Fisher's exact test for COG categories in pan, core genome analysis
I have come across several papers that do a Fisher's exact test to show over represented genes in COG categories specifically in a pan-core-accessory-species specific genome analysis. The test is done to show over represented genes in the core, accessory and species specific datasets as compared to the pan genome. Here is one such paper - Vibrio cincinnatiensis - comparative genome analysis

Usually, the way the pan and core genome analysis is represented in papers is by showing the proportion of genes in each of the COG categories in both. For example, the proportion of C (Energy production and conversion) could be 6.07 in the pan but 6.8 in the core genome. The Fisher's exact test is used to show that this particular category is over represented in the core genome as compared to the pan.

I already have the csv files from eggNOG from which I filtered out the genes in major COG categories for the pan, core, accessory and species-specific genomes into a single Excel sheet for 8 closely related genomes. Could someone tell me how to perform this test on this dataset?

My null hypothesis is that the proportion of genes in the pan, core, accessory and species specific genomes is the same. But, obviously this analysis is done to show how the proportions vary in some COG categories between the above mentioned 4 genome datasets.

I did try using python and R as well, but the test requires 2 rows and columns and here I am just comparing two columns (pan, core) with respect to 1 row - the 'C' category. It seems impossible to apply the test to this problem and I am sure I am missing something and that's why I posted it here.

I tried doing the Fisher's test on the proportions I mentioned previously (6.07 and 6.8) but I am getting an error : 'Fisher's test requires two rows and two columns'. Should I compare the absolute values, like '3745' genes in the pan and '4195' genes in the core?

# Answer
Sorry I've got my head around your question. For Fisher's Exact you need 2 "independent" observations and for SciPy its a 2x2 table.

So in your case ...

| Genome | Species A | Species B |
| --- | --- | --- |
| Pan | 6.07 | 6.10 |
| Core | 6.08 | 6.12 |

If you think about it, thats a better way of doing it because thats comparing pan vs core between species. A single observation isn't sufficient. To really do this you need 2 x 10 type contingency tables but SciPy will not do that, but you can do it in RStudio. For example using 10 species rather than just 2 ... and thats a real test of pan/core. My advice is to have your test species and in the other column have your control species ... like E. coli that sort of thing.

The SciPy solution is here. Place your data in a list of lists here its `mylist`. Next place into a numpy array and place that in a scipy function, like this...

```
import numpy as np
from scipy.stats import fisher_exact
mylist = [[6.07, 6.08], [6.10, 6.12]]
mytable = np.array(mylist)
res = fisher_exact(mytable, alternative='two-sided')
print (res.pvalue)

```

It should print the p value where \> 0.05 is not significant.

To start doing 2 x \>2 contingency tables in Python r2py is needed and from my experience (albeit about 5 years ago) I wouldn't advise it. r2py is a disaster when I tried to use it and just did the calculations in RStudio.
> 1 votes

---
Tags: statistics, eggnog, pangenome
---


title: Gene ensemble vs Hugo id

text:
# Question
Title: Gene ensemble vs Hugo id
Hi

I have been recently struggling with the gene annotations and changing gene ID to ensemble ID. While I was doing that there were some genes for the same species that had the same gene ID but different ensemble ID. Can someone let me know how this happens? for example, for these two ensemble IDs ENSCAFG00845015896, and ENSCAFG00000017519 which are both for canines there is the same symbol of "TNMD".

Thanks!

# Answer
> 2 votes
If you search for those gene IDs in Ensembl, you will see that ENSCAFG00845015896 is in the primary assembly, but ENSCAFG00000017519 is deprecated, and apparently has never been in the primary assembly but originates from the genome annotation of one individual dog. I did not align the protein sequences, but on a quick look they seem identical, which is why both are connected to the same gene symbol. If your identifiers are not all from the same assembly (and sometimes even if they are), you may see several ids for the same gene symbol.

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


# Question
Title: What is it about gene names starting with "LOC"?
I was struggling to use AnnotationDbi to change my ensemble ID to gene name, for datasets of three different species (human, canine, mouse). Among gene names for all three species there are genes with names starting with "LOC" like ENSG00000228037: LOC100996583, ENSCAFG00845012021: LOC613013 and ENSMUSG00000095523 LOC100038995.

What type of genes are they? and is the naming unique through species, I mean if I have LOC100996583 in humans, is it possible to have it in canines or mice? if yes, are they doing the same thing?

Thanks in advance!

# Answer
These are "genes" that don't have an official name. As explained in https://www.ncbi.nlm.nih.gov/books/NBK3840/#genefaq.Nomenclature:

> **Symbols beginning with LOC.** When a published symbol is not available, and orthologs have not yet been determined, Gene will provide a symbol that is constructed as 'LOC' + the GeneID. This is not retained when a replacement symbol has been identified, although queries by the LOC term are still supported. In other words, a record with the symbol LOC12345 is equivalent to GeneID = 12345. So if the symbol changes, the record can still be retrieved on the web using LOC12345 as a query, or from any file using GeneID = 12345.

In other words, when a gene has no known function or homologs in other species, when all we know about it is that it is a *locus* that seems to be actively transcribed, then it gets the name of `LOC` plus a numerical gene ID. So these genes, the LOC*NNNN* are essentially lesser known and studied genes. Some of them may graduate to a full gene name, others might not.

As for uniqueness, no absolutely not. GenBank gene IDs are unique, they don't imply homology. Only one gene can have a given ID. For example, gene ID 1 corresponds to human A1BG, alpha-1-B glycoprotein. The mouse homolog of human A1BG, mouse A1bg, alpha-1-B glycoprotein, has the gene ID 117586. The numerical value of the ID tells you nothing at all. It is just a unique identifier for a specific gene in a specific species.

Finally, a more general note, homology is rarely as straightforward as your question suggests you might think it is. Even actually homologous genes often have different functions with specific functions having been gained or lost in a given lineage. So it is very hard to know and never safe to assume that two homologous genes in different species "do the same thing".
> 3 votes

---
Tags: bioconductor, gene, genomics, human-genome
---


# Question
Title: What does the absence of a variant in a VCF file mean?
I have a `individual.snp.vcf.gz` file of an individual genome and the `referencegenome.snp.vcf.gz` file of the reference genome.

When I run the following code on the individual genome

```
gunzip -c individual.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

```

I get this:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  Type
chr1    10067   .   T   A   0   LowGQX  SNP
chr1    10070   .   C   A   0   LowGQX  SNP
chr1    10073   .   T   C   0   LowGQX  SNP
chr1    10105   .   A   C   0   LowGQX  SNP
chr1    10108   rs62651026  C   T   0   LowGQX  SNP
chr1    10157   .   T   C   0   LowGQX  SNP
chr1    10177   rs201752861 A   C   17  LowGQX  SNP
chr1    10180   rs201694901 T   C   0   LowGQX  SNP
chr1    10181   rs1246412344    A   T   0   LowGQX  SNP
chr1    10250   rs199706086 A   C   0   LowGQX  SNP
chr1    10257   rs111200574 A   C   0   LowGQX  SNP
chr1    10285   rs866375379 T   C   67  LowGQX  SNP
chr1    10321   rs1002315756    C   T   0   LowGQX  SNP
chr1    10327   rs112750067 T   C   36  LowGQX  SNP
chr1    10332   rs1175748383    C   A   0   LowGQX  SNP

```

When I run the equivalent code on the reference genome

```
gunzip -c referencegenome.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | head -n100 | awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

```

I get this:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  Type
1   10001   rs1570391677    T   A   .   .   SNP
1   10002   rs1570391692    A   C   .   .   SNP
1   10003   rs1570391694    A   C   .   .   SNP
1   10008   rs1570391698    A   G   .   .   SNP
1   10009   rs1570391702    A   G   .   .   SNP
1   10015   rs1570391706    A   G   .   .   SNP
1   10020   rs1570391708    A   C   .   .   SNP
1   10021   rs1570391710    A   G   .   .   SNP
1   10026   rs1570391712    A   C   .   .   SNP
1   10027   rs1570391716    A   C,G .   .   SNP
1   10032   rs1570391720    A   C   .   .   SNP
1   10033   rs1570391722    A   G   .   .   SNP
1   10039   rs978760828 A   C   .   .   SNP
1   10043   rs1008829651    T   A   .   .   SNP
1   10045   rs1570391729    A   C,G .   .   SNP
1   10051   rs1052373574    A   C,G .   .   SNP
1   10055   rs892501864 T   A   .   .   SNP
1   10056   rs1570391738    A   C   .   .   SNP

```

Now I want to get the genotype of the *individual genome* at all the sites listed in the *reference genome* table, so my question is: **When I don't find a particular variant ID in the *individual genome* table (for example, **rs1570391716**), can I assume that the genotype of the individual is equal to the homozygous *REF* (**AA**, for **rs1570391716**) in the *reference genome* table?** How certain can I be (assuming a low sequence reading error)?

# Answer
The absence of a variant means the variant was not found. However, that does not mean you can assume the position in the sequenced individual was homozygous for what is in the reference. Perhaps the position was simply not sequenced. Or the coverage was too low. So your first step will have to be to focus only on the regions covered by whatever sequencing assay was used to sequence your data, and only those regions with sufficient coverage for a variant to be called.

There is no magic number here, the threshold will depend on how the data were generated. For instance, a minimum coverage of 8 is a reasonable lower bound, but these days, we often have coverage values of 60x or more, and in the case of amplicon assays, we can have coverage in the thousands. So the threshold is something you will need to decide on.

Next, you really need to take the quality and filters into account. You haven't shown us the full VCF line, but my bet is that all the variants you show in your individual file are most likely not genotyped at all. These are not actually present. They all have a score of 0 and are marked with `LowGQX` which means:

> The genotyping quality locus (GQX) is less than 30 or not present.

And GQX is:

> The minimum of genotype quality assuming variant position and genotype quality assuming nonvariant position.

Variants with a QUAL score of 0 and low genotype quality are almost certainly not variants at all and you really shouldn't consider them.

So, since you apparently don't have the raw data, neither fastq nor bam, that means you cannot easily filter by coverage in the various positions. However, at the very least, you need to find out what sequencing assay and method was used to sequence your sample. Then, you can get a bed file with the regions of interest, and then limit your analysis to variants in those regions.

Also, I urge you to have a look at tools like `bedtools` and `bcftools` that will help you filter your data.

Finally, while your command will work, it is a bit inefficient, and you could simplify

```
gunzip -c individual.snp.vcf.gz | grep -v '^##' |
  awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' |
    awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

```

to

```
zgrep -v '^##'  individual.snp.vcf.gz |
  awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7, NR==1 ? "Type" : "SNP"}'

```
> 3 votes

---
Tags: vcf, genomics, genotyping
---


# Question
Title: Dealing with X Residues in FASTA files prior to folding analysis
We are working with a large dataset of proteins for folding/docking tests across various tools (AlphaFold, ESM2, RoseTTAFold). In some of the FASTA files, there is an X for a non-standard residue. An example of this is in the 6SEN\_1 FASTA file below. Is there a standard method of dealing with these non-standard residues prior to analysis to minimize impact on the folding/docking results?

```
>6SEN_1|Chains A, B|Transcriptional enhancer factor TEF-3|Homo sapiens (9606)
GPRSVASSKLWMLEFSAFLEQQQDPDTYNKHLFVHIGQSSPSYSDPYLEAVDIRQIYDKFPEKKGGLKDLFERGPSNAFFLVKFWADLNTNIEDEGSSFYGVSSQYESPENMIITCSTKVCSFGKQVVE**X**VETEYARYENGHYSYRIHRSPLCEYMINFIHKLKHLPEKYMMNSVLENFTILQVVTNRDTQETLLCIAYVFEVSASEHGAQHHIYRLVKE

```

# Answer
> 1 votes
Similar to @pippo1980 comment: start with a related standard amino acid, fold it however you are doing that, and then mutate to the nonstandard in the folded structure. You can do this with a toolkit like pdbfixer for standard AAs, you might need to do something like this to directly add the atoms for nonstandards, using an example pdb structure of the nonstandard AA.

You would then need to follow that with a round of energy minimization in the structure to resolve any large problems in the structure, using a toolkit like OpenMM.

If `X` truly means that you have no idea what's there, then... that's harder. In this specific case lysine may be a best bet, but it may be hard to generalize between e.g. methyllysine, phosphotyrosine, phosphoserine, etc. They're pretty different structurally, so having a good idea of what family of residue you would expect would be helpful.

The safest thing to do would be to go back to the transcript and in vitro translate to get an idea of what that X is. Though that's getting a little complicated compared to just replacing it with a leucine or whatever.

---
Tags: fasta
---


# Question
Title: Where can you download functional domain labels of the human genome?
I'm visualizing where mutations occur in a gene, and would like to add functional domains to the plots.

I have some trouble locating sources to download data on functional domains.

The dataset should contain:

* Functional domain name
* Start position in gene
* End position in gene
* Gene symbol name

Any tips kindly appreciated!

# Answer
> 0 votes
I would suggest looking into GFF3 annotations of the human genome.

If what you want is *protein domain annotations*, you may be more interested in the UniProt database which does seem to include genome annotation tracks for different protein functional regions.

I am unclear if uniprot directly has the gene names but this would be trivial to generate with a `bedtools intersect` operation.

---
Tags: genome, database, data-download, public-dataset, functional-annotation
---


# Question
Title: What does this BLASTN mean: Error pre-fetching sequence data?
I am trying to run BLASTN on a 14GB Pseudomonas database using the command: blastn -db ~/pseudomonas-db/pseudomonas -query ~/desktop/inputfiles/poz176.fasta -out poz176-out -num\_threads 8

I get the following error:

```
Error: NCBI C++ Exception:
    T0 "c++/include/corelib/ncbidiag.hpp", line 980: Critical: (CCoreException::eNullPtr) ncbi::CObject::ThrowNullPointerException() - Attempt to access NULL pointer.
    T0 "c++/include/corelib/ncbidiag.hpp", line 789: Critical: (CCoreException::eNullPtr) ncbi::BlastFormatter_PreFetchSequenceData() - Error pre-fetching sequence data
     Stack trace:
      /usr/lib/ncbi-blast+/libxncbi.so ???:0 ncbi::CObject::ThrowNullPointerException() offset=0xA6 addr=0x7fa1ea902f76
      /usr/lib/ncbi-blast+/libxobjmgr.so ???:0  offset=0x209609 addr=0x7fa1eb653609
      /usr/lib/ncbi-blast+/libxobjmgr.so ???:0 ncbi::objects::CBioseq_Handle::IsSetInst_Length() const offset=0x9 addr=0x7fa1eb653a09
      /usr/lib/ncbi-blast+/libxobjmgr.so ???:0 ncbi::objects::CBioseq_Handle::GetBioseqLength() const offset=0xE addr=0x7fa1eb653dce
      /usr/lib/ncbi-blast+/libxblast.so ???:0 ncbi::blast::LoadSequencesToScope(std::vector<ncbi::objects::CSeq_id_Handle, std::allocator<ncbi::objects::CSeq_id_Handle> >&, std::vector<ncbi::CRange<unsigned int>, std::allocator<ncbi::CRange<unsigned int> > >&, ncbi::CRef<ncbi::objects::CScope, ncbi::CObjectCounterLocker>&) offset=0x329 addr=0x7fa1ebd90e79
      blastn ???:0 ncbi::BlastFormatter_PreFetchSequenceData(ncbi::blast::CSearchResultSet const&, ncbi::CRef<ncbi::objects::CScope, ncbi::CObjectCounterLocker>, ncbi::blast::CFormattingArgs::EOutputFormat) offset=0x5FE addr=0x555af7fb31ee
      blastn ???:0 CBlastnApp::Run() offset=0xFDC addr=0x555af7fa99ac
      /usr/lib/ncbi-blast+/libxncbi.so ???:0 ncbi::CNcbiApplicationAPI::x_TryMain(ncbi::EAppDiagStream, char const*, int*, bool*) offset=0x193 addr=0x7fa1ea84edd3
      /usr/lib/ncbi-blast+/libxncbi.so ???:0 ncbi::CNcbiApplicationAPI::AppMain(int, char const* const*, char const* const*, ncbi::EAppDiagStream, char const*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) offset=0x61D addr=0x7fa1ea85032d
      blastn ???:0 main offset=0x6B addr=0x555af7fa7e5b
      /lib/x86_64-linux-gnu/libc.so.6 ???:0 __libc_start_main offset=0xEA addr=0x7fa1ea38ed0a
      blastn ???:0 _start offset=0x2A addr=0x555af7fa876a

```

I tried changing memory allocation but that didn't seem to do anything. I'm not sure what this error is telling me. How do I get around this?

(Had posted this on stackoverflow and was asked to post here)

# Answer
Based on other posts elsewhere it seems likely that there is something funky about your database. How did you generate it?

Could you try to re-generate the db without the `-parse-seqids` flag?
> 1 votes

---
Tags: blast, ncbi
---


# Question
Title: Variants from multiple tools: normalization before or after annotation with VEP?
We want to compare some variant calling tools and their calls on Whole Exome Sequencing data. We will have to normalize the variants that are called (.vcf format) before comparison using for instance vt normalize or hap.py. We are comparing with a gold standard too and use hap.py for this now to get some benchmark information.

Edit: with normalization I mean representing the variants in such a way, we are able to compare them: to ensure that the representation of a variant is “parsimonious” and “left-aligned.”

I was just wondering, and I cannot find any real guidelines (yet), when we annotate using VEP, do we have to normalize before or after annotation to keep the annotation correct as well. I know VEP has some normalization functionality too see here. So might we even not need an additional normalization tool before and just let VEP do it for us?

I am just wondering about the guidelines and what what others do. Thanks in advance,

# Answer
> 1 votes
If you're referring to left alignment + parsimonious representation, do it before any annotation. Most if not all annotation happens by comparing CHROM, POS and optionally REF and ALT and since vt normalize affects POS, you should do it before annotation.

I'm not sure how smart VEP is, but most annotators are simple lookup tools at least for third party databases, so they're not built to normalize while looking up.

---
Tags: variant-calling, vep
---


# Question
Title: Seeking Guidance on Designing DNA Sequences with Specific Criteria
Dear Bioinformatics Community,

I hope you're all doing well. I am currently working on a project that requires designing some 10 nucleotide single-stranded DNAs, (preferably), with specific criteria.

The criteria for the DNA sequences are as follows:

* An equal mix of nucleotides (A, C, G, and T)
* Thermodynamically stable sequences
* Sequences that do not self-hybridize or hybridize with each other

I have been exploring the NUPACK software suite to achieve these goals, but I am finding it challenging to use the tools effectively for this specific task. I understand that NUPACK is primarily designed for generating sequences based on target secondary structures, which may not directly align with my requirements.

I would greatly appreciate any advice, suggestions, or alternative methods you can provide to help me design these DNA sequences. Are there any other tools or resources that might be better suited for this purpose?

Thank you in advance for your insights and expertise. Your guidance will be invaluable as I work to complete this project.

Best regards,

# Answer
It is not difficult to simulate every 10mer DNA oligo, here is some Python code that will do that:

(there is probably a more elegant way to do this)

```
from itertools import product

nts = ["A", "T", "C", "G"]
# create a tuple with 10 copies of nts for the cartesian product
ten_tup = (nts,)*10

decamers = []
# create and iterate over a stringified cartesian product,
# which consists of all 10mers.
for a,b,c,d,e,f,g,h,i,j in product(*ten_tup):
  decamers.append(
    a+b+c+d+e+f+g+h+i+j
  )

```

At this point, you have some choices to make:

* how many 10mers do you want?
* what exactly is your criterion for GC content?
* what exactly is your criterion for thermodynamic interactions?

I'd suggest either primer3 or viennarna for evaluating thermodynamics. You can run each from the command line or from a web server, but they do have Python APIs, which I use below. ViennaRNA is definitely more trustworthy for sequences as short as 10, so I'll use it here. It will be a lot slower though as it's not cheating by using the nearest-neighbor approximations.

I have made up some numbers for the thresholds below, replace with your own.

```
from random import shuffle
GC_MIN = 4
GC_MAX = 6
DDG = -5
NUM_WANTED = 10

def test_gc(seq, gc_min=GC_MIN, gc_max=GC_MAX):
  gc_count = seq.count("G") + seq.count("C")
  if gc_count >= gc_min and gc_count <= gc_max:
    return True
  else:
    return False

def test_thermo(seq, other_seqs, ddg=DDG):
  intra = RNA.fold_compound(seq).mfe()[1]
  print(intra)
  if other_seqs:
    # find the most favorable interaction with any other sequence
    inter = min(
      [RNA.fold_compound(seq+"&"+other_seq).mfe_dimer()[1] for other_seq in other_seqs]
    )
  else:
    inter = 0
  homo = RNA.fold_compound(seq+"&"+seq).mfe_dimer()[1]
  if intra > ddg and inter > ddg and homo > ddg:
    return True
  else:
    print("thermo reject!!")
    return False

chosen_seqs = []
# shuffle the decamers to ensure that you get dissimilar sequences
shuffle(decamers)
i = 0
for seq in decamers:
  i+=1
  print(i)
  if test_gc(seq) and test_thermo(seq, chosen_seqs):
    chosen_seqs.append(seq)
  if len(chosen_seqs) >= NUM_WANTED:
    break

# chosen_seqs has your decamers that pass the threshold

```

In principle, you could even combine these 2 steps but I think it's cleaner to keep them apart unless storing 4^10 ~= 1 million decamers in memory is an issue, which it shouldn't be with modern computers.

**Caveat emptor, this code may contain bugs and I haven't tested it other than making sure it runs.**
> 1 votes

---
Tags: computational-biochemistry, design, dna, thermodynamics
---


# Question
Title: How to normalise Bulk RNA-seq data to account for transcript length, coverage depth, library size and batch effects?
I'm currently working on a project where I'm comparing aggregate measurements (mean, median, etc.) of expression data (RNA-seq) from different groups of genes across various samples with different characteristics (tissue type, health status, etc.). Additionally, the raw counts were collected from several different labs using various techniques.

Since I am conducting between-gene measurements, the data should be normalised to account for differences in transcript length and coverage depth (TPM, RPKM, FPKM). However, I am also interested in comparisons across samples based on tissue type and other factors. Therefore, the data should also be normalised to account for library size (TMM, quantile, etc.), and, as the data were collected from multiple sources, it should be corrected for batch effects.

I have read through many papers but am unsure and confused about how to proceed with the normalisation procedure starting with the raw counts. Can I simply string the methods together, starting with batch effect correction, followed by library size normalisation, and then the within-sample normalisations?

I would appreciate any insights or suggestions on this. Thanks

# Answer
> 2 votes
All factors you mention bar the batch effects can be addressed by some sort of normalization (I use a modified version of UQpgQ2). Batch effects can be addressed using ComBat/SVA/ComBat-seq but my team's experience with them has been bad - they introduce artifacts and skew the data quite a bit especially if the batch effect is not pronounced.

Unfortunately, there is no way to do simple comparisons without huge caveats. You're better off using DESeq2 to perform statistically sound comparisons using the batch variable as a covariate.

# Answer
> 1 votes
Have a look through the DESeq2 manual. It goes into a lot of detail about how it works, and how you can use it to do all the things you're asking about:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

The standard workflow already has all these things strung together, and works fairly well for the most common differential expression tasks:

> # Quick start
>
> Here we show the most basic steps for a differential expression analysis. There are a variety of steps upstream of DESeq2 that result in the generation of counts or estimated counts for each sample, which we will discuss in the sections below. This code chunk assumes that you have a count matrix called cts and a table of sample information called coldata. The design indicates how to model the samples, here, that we want to measure the effect of the condition, controlling for batch differences. The two factor variables batch and condition should be columns of coldata.
>
> ```
> dds <- DESeqDataSetFromMatrix(countData = cts,
>                              colData = coldata,
>                              design= ~ batch + condition)
> dds <- DESeq(dds)
> resultsNames(dds) # lists the coefficients
> res <- results(dds, name="condition_trt_vs_untrt")
> # or to shrink log fold changes association with condition:
> res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")
>
> ```

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


# Question
Title: Can anyone share insights or provide the complete mRNA sequence for tdtomato-WPRE for scRNA seq analysis?
I'm currently immersed in research involving the Rosa26-tdtomato (Ai9) mouse line. We performed scRNA seq and analysis focusing on tdtomato-expressing cells but faced a challenge with a lower-than-expected number of positive cells.

Our alignment process, using the tdtomato sequence (Genbank accession: AY678268), seems insufficient. Fellow researchers suggest success in similar cases by incorporating the full mRNA sequence, including the 3´UTR, along with tdtomato-WPRE mRNA.

I searched everywhere for the full mRNA sequence but couldn't find it. Can anyone share insights or provide the complete mRNA sequence for tdtomato-WPRE? This would greatly enhance our alignment process and improve the accuracy of tdtomato positive cell detection.

Appreciate your contributions to advancing our research!

# Answer
> 2 votes
I'd suggest looking at the mouse itself, which appears to be here on the JAX website.

This link includes RNA-seq experiments, so presumably you should be able to compare your pipeline and data to previous work to see whether there's anything weird there.

Looking at the original paper, you can see the structure of the transgene which might be enough info to work from, there is also information in the supplement. You should in principle be able to figure out the sequence of each of those components and include them.

If you are really having issues, you could do a transcriptome assembly to see whether you can recover the locus.

You can always write to the senior author asking about this line, or try to figure out why your transcript may truly not be expressed.

# Answer
> 0 votes
The Genbank sequence you've specified appears to be for RFP, not tdTomato.

Have you tried one of the other tdTomato sequences, for example a lentiviral sequence, or a sequence that is called tdTomato?

https://www.ncbi.nlm.nih.gov/nuccore/KJ081243.1

---
Tags: scrnaseq, seurat, tdtomato
---


# Question
Title: WGCNA: Extract genes within the different modules to create a dataframe with two columns
My inquiry pertains to WGCNA. After identifying the modules, I would like to extract genes within the different modules to create a dataframe with two columns: `Gene_a` and `Gene_b`. The values in these columns should represent co-expressed genes. An illustrative example is as follows:

In this example, the gene `XX62730x` is co-expressed with `XCGR1Y`. And `YYYYY` is co-expressed with `ZZZZZZ`.

However, I am facing difficulties as the official WGCNA page is no longer accessible.

I would greatly appreciate any help you can provide. Thank you in advance!

Best,
Nana

# Answer
You can get this information, more or less, using the functions exportNetworkToVisANT and exportNetworkToCytoscape. You will find some extra columns that you can ignore or delete from the resulting table. We are still working with UCLA to restore the WGCNA page including the tutorials, but the tutorials can temporarily be accessed using this dropbox link.

To use either one of the two functions, you will need the network matrix for which we usually use TOM. You will probably want to set a threshold for connections to be included, otherwise you will get a fully connected network and the resulting table will (usually) be huge.

Hope this helps!
> 1 votes

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


# Question
Title: Discrepancy in Depth of Coverage Estimation between Pileup and samtools depth
I'm encountering a discrepancy in the depth of coverage assessment for specific coordinates when comparing results obtained from pileup files and samtools depth.

I've processed pileup files and extracted coverage information for a given coordinate, resulting in a coverage depth of **8028**. However, when using samtools depth on the same coordinate, the coverage depth is reported as **88565**.

My analysis involves the following steps:

Alignment of paired-end reads using bwa mem. Generation of pileup files and extraction of coverage information. Utilizing samtools depth to calculate coverage depth for the same coordinates. Visual inspection in IGV for validation, where I observe that the number of reads aligning to the region matches the count obtained from samtools depth. I've noticed that in this region, there are no structural variations. The versions of samtools I'm using are samtools 1.19.2.

Could you please help me understand the potential reasons for such big discrepancies? What parameters should I review, and could the differences be attributed to any specific factors?

Thank you for your assistance!

# Answer
This is purely a guess based on what you write given that you are not showing the exact command for samtools etc.

I suspect that you are allowing sequences with MAPQ==0 (supplementary alignments) into your depth calculation. A brief google search tells me that this issue has been previously observed with `samtools depth`.

In principle, you should be able to see this directly in IGV, by controlling how supplementary alignments are displayed.

This is something that you can directly control with the `-Q` flag, to restrict it to uniquely mapping reads (where you set the flag to e.g. `-Q 1` or higher).

I'd suggest trying to play with `samtools depth` parameters and seeing if you can bring it more into alignment.
> 1 votes

---
Tags: samtools, coverage, mpileup
---


# Question
Title: Fastest way to convert GTF or GFF file to a BED file?
Given a GTF/GFF3 annotation file, what is the fastest way to convert it into a BED file?

# Answer
# gxf2bed

The fastest G{F,T}F-to-BED converter around the block!

translates

```
chr27 gxf2bed gene 17266470 17285418 . + . gene_id "ENSG00000151743";
chr27 gxf2bed transcript 17266470 17281218 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8";
chr27 gxf2bed exon 17266470 17266572 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8"; exon_number "1"; exon_id "ENST00000541931.8.1";

```

into

```
chr27 17266469 17281218 ENST00000541931.8 1000 + 17266469 17281218 0,0,200 2 103,74, 0,14675,

```

before your eyes blink!

Converts

* *Homo sapiens* GRCh38 GENCODE 44 (252,835 transcripts) in 3.34 seconds.
* *Mus musculus* GRCm39 GENCODE 44 (149,547 transcritps) in 2.11 seconds.
* *Canis lupus familiaris* ROS\_Cfam\_1.0 Ensembl 110 (55,335 transcripts) in 0.99 seconds.
* *Gallus galus* bGalGal1 Ensembl 110 (72,689 transcripts) in 1.14 seconds.

## Usage

```
Usage: gxf2bed[EXE] --input/-i <GTF/GFF> --output/-o <BED>

Arguments:
    --input/-i <GTF/GFF>: a .gtf/.gff file
    --output/-o <BED>: path to output .bed file

Options:
    --help: print help
    --version: print version
    --threads/-t: number of threads (default: max ncpus)

```

## Installation

to install gxf2bed on your system follow this steps:

1. get rust: `curl https://sh.rustup.rs -sSf | sh` on unix, or go here for other options
2. run `cargo install gxf2bed` (make sure `~/.cargo/bin` is in your `$PATH` before running it)
3. use `gxf2bed` with the required arguments
4. enjoy!

## Build

to build gxf2bed from this repo, do:

1. get rust (as described above)
2. run `git clone https://github.com/alejandrogzi/gxf2bed.git && cd gxf2bed`
3. run `cargo run --release -- -i <GTF/GFF> -o <BED>`

## Container image

to build the development container image:

1. run `git clone https://github.com/alejandrogzi/gxf2bed.git && cd gxf2bed`
2. initialize docker with `start docker` or `systemctl start docker`
3. build the image `docker image build --tag gxf2bed .`
4. run `docker run --rm -v "[dir_where_your_gtf_is]:/dir" gxf2bed -i /dir/<GTF/GFF> -o /dir/<BED>`

## Conda

to use gxf2bed through Conda just:

1. `conda install gxf2bed -c bioconda` or `conda create -n gxf2bed -c bioconda gxf2bed`

## Benchmark + FAQ

If you google "gtf to bed" or "gff to bed" you'll find some posts about recommended tools. Part of these tools are already deprecated or too difficult to run (even in Linux), due to poor maintenance. This project was conceived to provide an easy way to convert GTF/GFF files to BED structures and finish a set of high-performance converters between gene model formats in Rust (bed2gtf, bed2gff and now gxf2bed).

To test the efficiency of gxf2bed, I took 4 GTF/GFF-to-BED converters and ran them with the current GRCh38 GENCODE annotation (v.44), that has ~250,000 transcripts \[the values showed here are the mean of 5 consecutive runs\]. Briefly, the competitors are:

* UCSC's utils: `gtfToGenePred | genePredToBed` & `gffToGenePred | genePredToBed` (taken from here)
* Signal & Brown: `gtf2bed.py` & `gff32gtf.py && gtf2bed.py` (taken from here)
* ea-utils: `gtf2bed.pl` & `gff2gtf.pl | gtf2bed.pl` (taken from here)
* BEDOPS: `gtf2bed` & `gtf2bed` (taken from here)

Besides the easy way to make `gxf2bed` run and the GFF-GTF channel that allows centralize both in a single converter, this tool showed a significant difference in computation time against the other tools in each one of the two formats. Even using a combination of different number of threads, times were practically the same 3.5s +/- 0.4s and the differences were maintained. Is important to note that GFF3 files took a lot more time than its GTF counterpart. This could be due to the need of chaining two different programs 1) GFF-to-GTF and 2) GTF-to-BED.

Taken together, `gxf2bed` offers an easier and faster way to convert GTF/GFF files to BED files. This tool could save you at least x2-3 times the time you used to spend using other tools to convert GTF files and x5 times if you want to convert GFF3 files!

---

Disclaimer: I am the author of the gxf2bed tool.
> 3 votes

---
Tags: bed, format-conversion, gtf, gff
---


# Question
Title: split fastq file containing a sequence block at different locations
I have some fastq files (obtained from nanopore sequencing) that contain reads that can be of either of these 5 forms:

* a known CDS with 3'UTR:

```
CDS----------------------- (seq1: original sequence)

```
* the same cds with a block of 50bp represented by `|||||` (that substitutes the 3'UTR at 4 different locations)

```
CDS|||||------------------ (sub1: substitution 1)
CDS-----|||||------------- (sub2: substitution 2)
CDS----------|||||-------- (sub3: substitution 3)
CDS---------------|||||--- (sub4: substitution 4)

```

I need to split my fastq files in 5 smaller files clustering the different types of sequences (seq1, sub1, sub2, sub3 and sub4). The problem is that those sequences all share some segments and I don't know which method I should apply. I tried finding the best match by building a BLAST database with the 5 sequences but as it is local alignments, I get HSPs for each sequences and cannot discriminate between them.

Does anyone has an idea how I could achieve that?

Any help would be amazing.

**edit**

Here are 5 sequences as I expect them (i.e., fasta reference)

```
>seq1 (this one does not have the 50bp block)
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub1
_GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub2
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACT _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub3
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACT _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub4
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTA _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ GATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA

```

**examples (50bp block if present is in lowercase here)**

* this read should be associated with sub3 because of the flanking sequences (...CTACAACT in 5' and TTGTCTCTTTT... in 3' of the 50 bp block)

```
@e213fdd5-bbfb-4dfe-9b03-29eea256da0d
TGTGTACTTCGTTCAGTTACGTATTACTAGTTATTGAGTGTCTTTGTGTTTCTGTTGGTGCGTCTTCGCACAAGGCTAATCTTACTCTCCTCTCCAATGTACTGCTTCTATGTTCATCTCAGCCGAATTCAATAAGGAGAAGAACTTTCACTAGGTTTCCCCATTCTGTTCGATTAATTGGTGATGTTAGTAAATTTCAATTTTCTGTCGGTGGAAGGTGAAGGTGATTGCTGCTGGCTTTTCATTGTTTACTTTTGGCTTCTGTTTTGTAGTAACGCATCACTACTTTCTCTTATAGTGTTCAATGCTTTCCAAGTCTGGAATCATGAAGCAGCGGCTTCTTCAAGAGCGCCATGCCTGAGGGATACGTGCAGGAGAGGACCATCTTCTTCAAGGACGACGGGAACTACAAGACGTGCTGAAGTCAAGTTTGGGGAGACACCCTCGTCAATGAGATCGAGCTTAAGGGAATCGATTTCAAGGAGGACGGAAACATCCTCGGCCACAAGTTGGAATACAACTACAACTCCCACAACGTATACATCATAGCCGACAAAACCCAAAAGACGGCATCAAGAAACCAACTTCAGACCTTACCATCAGGCGGCGGCCATCACACAGTCATTGTCAACAAATGTAATTGGCGGTGGCTATTCTTTTTACCCAGACAACCATTACCTTGTCCACACAGTCACACTGTGAAAAGATCCCAACTTTGAAAGAGAGGCCACTGGTCTTCTTGAGTTTGTAACAACTGCTGGGATTACACATGGCATGGATGAACTATACAAACATGACAGACTCTAAAGCCTGCAGTTCCTGAAGAGATTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTgtaagtgagggagtctcggactcaattagaggctcttctttcacaTTGTCTCTTTTGGTTTATTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTTTCCTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATGCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCACCATCAATGAAGATAGAGCGACAGGCAAGTCACAAAAACACCGACAACTTTCTTGTCACGGTAGGCGATAGCAATACGTAACTG
+
$&&&%234441124--585545;>B99968650-++*()'&'(;D540488?=888,++*+*)*$$$##&%$$$%)*****('%$$'1/+++-,('))'%$$%$$$&)&%%%&(*'&),/'''&&'06872..*&%)),,3/49:;<?A773)('()/066+'&&)('(+&%$&''$$&&()*--888)('()'')')'''++*+)&')&$$$%()'')*''&&))&'&$&(&$#$%&%&&&%%$%&..&%%%%&&%$%&'&$$'-)(%%%$&(*,+****&'''('+17<{..+*('((*+)*&&&'**.*%$&%&&%$$$%&&'&&(+'(&&&%%&&($$%&&'+55536988;;74632+(*)''%6<=(((&*''(68998===?ABBBB=B?<;;96667>>?>>?@?A=-+.7556<@>A<0//002(()-989::::;<<>>;:5++++->;::;<?F:7<()988;;<?>9;9)36:5-963/002B@A@@?554543,,9<=<@?@>777461166.../B88<;<@=<:888;;>42(''%%(%%+*('&%%&&**)(()*+020*'('&$$$'('''&%'((%(&$##$&)*()&%%%$$#$%'&(&'%$$%%&'%%%$%%%$&('$##$&&&&$$$%$%$$%&%%&040(((&'.023266,,,++0)++-.-*,,0+&&'$&''''%%'&-+,0((+))*+*(%&)+00.3)''%'&'''(*(%%)++++,--+.....0'(()7:<=;=<<<<<0.<;9:;<@;;;:<<=;2+)++,1212611)())*:1++*'&&%((&&056/,454))*''''(3346;<;;;;;000.0/043546??=;<*)+64335663455665333788:78779<>9755446::<><198662*+7899722248867>?8++4))+6779<<BA><9;<;:;;{<?BB-,,-13444/154///933{3==@B:9989::::::;;<<=:3335=EBB;:*)..'-:69?@CKKKDC<<?>:)()((012.,,+,55;=?>@>:98:78544344*)%$$$&)))))))))(((((''''&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')+.1866677978:;;:::;;;99::988889653456:::9889;::8899.,+0668:;:::8436735//:632/-,+**+013459;9332/+($

```
* this read should be associated with sub4 because of the flanking sequences (...TTTTACTCCTA in 5' even if a C is missing, and GATCTTGTCCT... in 3')

```
@107dfd6d-2101-4fd8-aeb5-2c432e176d72
ATGTACTTCGTTCAGTTACGTATTGCTATCGCCTACCGTGACAAGAAAGTTGTCGGTGTCTTTGTGTTTCTGTTGGTGCTGATATTGCATGAAGACTAATCTTTTCTCTTTCTCATCTTTTCACTTCTCCTATCATTATCCTCGACCGAATTCAGTAAAGGAGAAGAACTTTTCACTGGGTTGTCCCAATTCTTGTTGAATTAGATGGTGATGTTAATGGGTACAAATTTTCTGTCAGTGGAAGAGGGTGAAGGTGATGCAACATACGGAAAACTTACCCTTAAATTTATTTGCACTACTGGAAAACTACCTGTTCCATGGCCAACACTTGTCACTACTTTCTCTTATGGTGTTCAATGCTTTCAAGATACCAGATCATATGAAGCGGCACGACTTCTTCAAGAGCGCCATGCCTGAGGGATACGTGCGGGAGAGGACCATCTTCTTCAAGGACGACGGGAACTACAAGACACGTGCTGAAGTCAAGTTTGAGGAGACACCCTCGTCAACAGGATCGAGCTTAAGGGAATCGATTTCAAGGAGGACGGAAACATCCTCGGCCACAAGTTGGAATACAACCGCAACTCCCACAACGTATACATCATGGCCGACAAGCAAAAGAACGGCATCAAAGCCAACTTCAAGACCCGCCACAACATCGAAGACGGCGGCGTGCAACTCGCTGATCATTATCAACAAAATACTCCAATTGGCGATGGCCTGTCCTTTTACCAGACAACCATTACCTGTCCACACAATCTGCCCTTTCGAAAGATCCCAACGAAAAATTTTTGCATGGTCCTTCTTGAGTTTGTAACGGCTGCTGGGATTACATGGCATGGATGAAGAACTGCAACACAAACACTTGACGAACTCTAAACCTGCCAGTTCCTGAAGAAATTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAAGTGACGATGGAACCTGTGTCACCAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCTAgtaagtgagggagtctcggactcaattagaggctcttctttcacaGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCACCATCAATGAAGATAGAGCGACAGGCAAGTCACAAAGACACCGACAACTTTCTTGTCACGGTAGGCGATAGCAATACGTAACT
+
-0,+))),:::<??===>=<<<?@>=<<=?;:889:<;<<=;;;;CD@?:9:A>9999=999DBAA@<BE=3220111<B@@AA@A?>>>=>>;:9999:9;<E9A<;?????C8788;==@?;;;==9989@A6677?=7788+**+368>B87::==@<99;<<<<<<??>@=:9/.))*698.*,/8A=<==>>=?BAB=::::;;;<?@>@B@@?:888:>?@CCCEDB>??>===::.+)()2157=7:99<B?>???@??=;::=DGC=<6?;84568>AB@CCBCBAA====AAACEI{AA=<::;<?@=<<>=>@AAA@>>>>?<?====A?==<=@>??@>?==??CA@AAA@AB=757:965(*.;<;<::?AAA@>=>::::<=?@BCC5336899=?>>?@A??42220;9:89.-***+77;97<;:;<<??BC{878<:9:<:2112<<>===@=<;;<>???877/,.-.-///0,++,)299:=;:;;=@@<;=7778?=::;;<;=<<C>AB<BBBCEDCDC?A>;;9:655:DAA@?>>>=?=><;<<@=>><=B=@=:870++,,;=;>><<==CBA@?>=>>=667<;;7778B>B?B<<==<99:::4:659;6642212:.-<=??=6556;??@@>::;:==;8887899==;2143235<<<?>==?CEBD?>>>ADB=720/.02334<00//06568>:888:711*)-776'''&06:<=><><866;{C@A?>8897=>EK558568552-179664)(()),,--,()))3:;:<BAAA>=>>>><<>>*))+/245999;>>=96;?B9888<<77-,,,,&&%%%&&'+6;A;2111--,,.:;><<ACCA@@?=9;;;>=>==@930../+)(())07>>8<><;:;9238<>=>998878:7700016545**1236434689<9<<;.--0-,,,&'():78862237@???><;;;;<<=@>@??@@===;<=>>;:8;)(67511223/.---.2001023:;21<?=**6223;;>?AACGA=73366:;<=:766;ACC=<;86000?8=@DHEC@@=;;<9999<<==9><==8889>A>=;30266777:972+***)))'&%%%%%$$$$$$######"""""""""""""""""""""""""""""""""""""""""#&')+-1122/...46688<<=0('%%&258899;92136:999;:::9:;;--2-.343668=<?@EF????==:7668;6667;><<;;>;;8851.*

```

# Answer
You can try seqkit locate with say start or end of CDS sequence. This will give you a tabulated output with locations of your CDS query in each read. Since you expect that the CDS sequence should be located at one of the five different locations inside your reads you can figure out to which one of the five a given read belongs.

You may also use clumpify.sh from BBMap to cluster reads prior to searching for the CDS locations. The idea is that maybe the five groups will cluster in fairly well separated blocks, but this is fairly speculative.

*edit*

Since the alignment of your sequences looks like this:

```
CLUSTAL 2.1 multiple sequence alignment

sub1            NNNNNNNNNNGTAAGTGAGGGAGTCTCGGACTCAAT---TAGAGGCTCTTCTTTCACANN
sub2            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTNNNNNNNNNNGTA
seq1            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGG-
sub3            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGC
sub4            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGG-
                          **   *** **     * **  *  *   *  *  **

sub1            NNNNNNNNTTGATCCAACGGCTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
sub2            AGTGAGGGAGTCTCGGACT-CAATTAGAGGCTCTTCTTTCAC---ANNNNNNNNNNTTTG
seq1            --------------------CTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
sub3            TAGAAGTGACGATGGAACCTGTGTCTACAACTNNNNNNNNNNGTAAGTGAGGGAGTCTCG
sub4            --------------------CTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
                                              *                          * *

sub1            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
sub2            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
seq1            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
sub3            GACTCAATTAGA--GGCTCTTCTTTCACANNNNNNNNNN--TTGTCTCTTT-----TGGT
sub4            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTANNNNNNNNNNGTAAGTGAG
                  *   ****    *  * **  **  *                            **

sub1            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub2            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
seq1            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub3            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub4            GGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACANNNNNNNNNNGATCTTGTCCTCTT
                  * * *                *** * * *   *          **************

sub1            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub2            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
seq1            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub3            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub4            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
                ************************************************************

sub1            ATAAA
sub2            ATAAA
seq1            ATAAA
sub3            ATAAA
sub4            ATAAA
                *****

```

instead of CDS end/start you have to select regions around your deletions with say 15-20bp flanks on each side. Put these in fasta and run locate.

**caveat** running with zero mismatches will miss some reads.

**edit 2**

The command:

```
zcat input.fq.gz | seqkit locate --threads 8 \
--max-mismatch 0 \
--pattern-file fasta_with_5_frags_extracted_from_algnt.fa > result.tsv

```

Some trial and error is required to select good fragments matching exactly one of the five groups plus increasing `max-mismatch` if due to sequencing errors a lot of reads were dropped. To speed up IO one can compress `result.tsv` on the fly, i.e. using `pigz` or `zstd`.

In the next step we need to create five separate files listing read names containing a given fragment:

```
grep frag_name_1 result.tsv | cut -f 1 > frag_name_1.read_names.txt

seqkit grep --pattern-file frag_name_1.read_names.txt input.fq.gz -o frag_name_1.reads.fq.gz

# repeat five times with different pattern file and changing the output

```
> 4 votes

---
Tags: sequence-alignment, ngs, clustering
---


title: How to classify TF motifs by family

text:
# Question
Title: How to classify TF motifs by family
I have a FIMO output, the input sequences were putative promoter sequences of several species. I want to graph the positions of these motifs along a horizontal graph, but I want to filter the data by clustering the motifs by family, e.g. all the zinc finger motifs being represented by one color, all the helix loop helix motifs being represented by one color, etc. How would I go about clustering the raw .tsv file?

# Answer
> 0 votes
Comment: I don't understand what you mean by "I want to graph the positions of these motifs along a horizontal graph"? Are you talking about a Manhatten plot? Or something else? And I'm not sure about "filter the data by clustering the motifs by family". What type of clustering do you want to perform?

Without further information/description it's difficult to understand what you're trying to do. If possible, please edit your question with an example of what your desired output. Also please indicate the language you're using for your visualisation, e.g. R, python, other.

---

Answer: I may have completely misunderstood, but perhaps this approach could be adapted to suit. Using the example FIMO output file from https://meme-suite.org/meme/doc/examples/fimo\_example\_output\_files/fimo.html?man\_type=web:

```
library(tidyverse)
library(qqman)
#>
#> For example usage please run: vignette('qqman')
#>
#> Citation appreciated but not required:
#> Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731
#>

fimo_df <- read_tsv("~/Desktop/fimo.tsv", comment = "#") %>%
  rename("SNP" = motif_id, "CHR" = sequence_name, "BP" = start,
         "P" = `p-value`, "zscore" = `q-value`) %>%
  mutate(CHR = case_when(str_detect(CHR, "X") ~ "23",
                         str_detect(CHR, "Y") ~ "24",
                         TRUE ~ CHR)) %>%
  mutate(CHR = parse_number(CHR)) %>%
  select(-c(stop, strand, matched_sequence))
#> Rows: 1494 Columns: 10
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr (5): motif_id, motif_alt_id, sequence_name, strand, matched_sequence
#> dbl (5): start, stop, score, p-value, q-value
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

manhattan(fimo_df, p = "zscore", logp = TRUE, ylab = "q-score",
          genomewideline = FALSE, suggestiveline = FALSE,
          main = "Manhattan plot of TF motifs",
          chrlabs = c(1:19, "X", "Y"))

```

If you want to colour the points by TF motif (or by 'cluster'), perhaps ggplot would be appropriate, e.g.

```
library(gghighlight)
axis_set <- fimo_df |>
  group_by(CHR) |>
  summarize(center = mean(BP, na.rm = TRUE))

ylim <- fimo_df |>
  filter(P == min(P)) |>
  mutate(ylim = abs(floor(log10(P))) + 2) |>
  pull(ylim)

ggplot(fimo_df, aes(x = BP, y = -log10(P),
                    color = as.factor(CHR),
                    size = -log10(P))) +
  geom_point(alpha = 0.75) +
  # scale_x_continuous(label = sort(axis_set$CHR), breaks = sort(axis_set$center)) +
  scale_y_continuous(expand = c(0, 0), limits = c(3.95, ylim)) +
  scale_color_manual(values = rep(
    c("#276FBF", "#183059"),
    unique(length(axis_set$CHR))
  )) +
  scale_size_continuous(range = c(0.5, 3)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
  ) +
  gghighlight(motif_alt_id == "TBP")
#> Warning: Tried to calculate with group_by(), but the calculation failed.
#> Falling back to ungrouped filter operation...

```

Created on 2024-02-05 with reprex v2.1.0

Does that help at all?

# Answer
> 0 votes
How do you know the "family" to which your motifs belong? If the source for the motifs that FIMO matched is, say, JASPAR Vertebrata CORE, then you might reasonably decide to use JASPAR's "structural class" as provided by their pre-computed "Annotation table" in JASPAR CORE Vertebrates clustering.

Assuming you are using JASPAR and you do make this decision...

Your "raw .tsv file" is in FIMO Results TSV Format. "Joining" it with a downloaded copy of the "Annotation table" serves to annotate it with additional columns, having the effect you seek of clustering the FIMO results on the "structural class" of each matched motif. Joining can be implemented following Running queries directly against CSV or JSON using the sqlite-utils command-line utility.

How you would then plot these features and encode the class as a color requires more detailed description what the plot is supposed to represent or look like. For instance, do you "want to graph the positions of these motifs along a horizontal graph" with coordinates taken with respect to the TSS of the transcript from which the promoter is taken?

---
Tags: sequence-analysis, motifs
---


# Question
Title: Finding homologous regions in multiple whole genomes
Right now I have six genomes that I want to compare and identify homologous regions in the genomes. I have run nucmer, show-coords and obtained the output files. An example is shown below with Genome 1 vs Genome 2. More files go through Genome 1 vs Genome 3 until all are compared.

I was trying to use bedtools merge after concatenating the multiple output files from show-coords command to find the homologous regions. However, I am not quite sure on the command to use. Any help is appreciated.

Edit: added commands. I just want to know what regions are homologous across the six genomes (99/100% IDY). I don't know what I am supposed to merge to find this out. I was reading and I believe bedtools should be but I don't know. I am usually a wetlab but I try to do a bit of bioinformatics.

```
 # Loop through all pairs of genomes
for i in {1..6}; do
    for j in $(seq $((i + 1)) 6); do
        genome1="genome${i}.fasta"
        genome2="genome${j}.fasta"

        output_prefix="output_${i}_${j}"

        # Run NUCmer
        nucmer -p "${output_prefix}" "$genome1" "$genome2"

        # Show coordinates
        show-coords -r -l -c -T "${output_prefix}.delta" > "${output_prefix}_coords.txt"
    done
done

```

After using the above I then concatenate the files with

```
cat output_*_coords.txt > combined_coords.txt

```

This gives me something that looks like this below.

```
[S1]    [E1]    [S2]    [E2]    [LEN 1] [LEN 2] [% IDY] [LEN R] [LEN Q] [COV R] [COV Q] [TAGS]
1       14347   1       14347   14347   14347   99.02   5064019 5004926 0.28    0.29    Genome1      Genome2
14474   32175   14350   32051   17702   17702   99.33   5064019 5004926 0.35    0.35    Genome1      Genome2
16541   22042   4698474 4692975 5502    5500    99.64   5064019 5004926 0.11    0.11    Genome1      Genome2
16548   22042   3728461 3722968 5495    5494    99.84   5064019 5004926 0.11    0.11    Genome1      Genome2
16583   22041   3904238 3898781 5459    5458    99.65   5064019 5004926 0.11    0.11    Genome1      Genome2
18727   22097   4426876 4423506 3371    3371    99.73   5064019 5004926 0.07    0.07    Genome1      Genome2

```

# Answer
> 2 votes
According to show cords the format you are getting contains columns which you can use to obtain `BED` format with relevant information.
You need following columns:

```
TAGS_genome1
[S1]    Start of the alignment region in the reference sequence.
[E1]    End of the alignment region in the reference sequence.
TAGS_genome2
[% IDY]

```

I would start with not merging all coordinates but getting bed files for each pairwise genome\_a\__vs\__genome\_b.bed etc. and start from there, taking a subset of 3 genomes first (say a, b,c). You can find the regions shared by all three `bedtools intersect`:

```
bedtools intersect -a a_vs_b.bed -b a_vs_c.bed

```

The issue is that you will have many to many relationships, since region\_1 from b can have several similar regions in b, and each one of these (== regions in b) also can be similar to a number of regions in a. Depending on what you are looking for and how dis-similar the genomes are you may end up with reporting repetitive elements or if you filter hard missing i.e. exons.

---
Tags: sequence-alignment, bedtools
---


# Question
Title: How to calculate distance between two regions from bed like dataframe in R?
I'm wondering how can I calculate the distance between two regions from the bed (`chr`, `start`, `end`) alike dataframe in R. Precisely, I want to subtract `end` of r from `start` of r+1 (r=row). Also, the calculation will repeat similarly for each chromosome.

```
head(G3_dat[1:3])

# A tibble: 6 × 3
  chrom   start     end
  <chr>   <int>   <int>
1 chr1  1330426 1330628
2 chr1  2017424 2017750
3 chr1  2017424 2017750
4 chr1  2018462 2018871
5 chr1  3899411 3899468
6 chr1  4653431 4653724

```

The expected outcome would be,

```
head(G3_dat[1:4])

    # A tibble: 6 × 4
      chrom   start     end     dis
      <chr>   <int>   <int>     <int>
    1 chr1  1330426 1330628
    2 chr1  2017424 2017750     686796
    3 chr1  2017424 2017750       -326
    4 chr1  2018462 2018871        712
    5 chr1  3899411 3899468    1880540
    6 chr1  4653431 4653724     753963

```

# Answer
There are two functions in the tidyverse suite (specifically `dplyr`) called `lead` and `lag` that allow calculations like this to be done:

```
library(tidyverse)
G3_dat %>%
   group_by(chrom) %>%
   arrange(start) %>%
   mutate(dis = start - lag(end))

```

I'm not sure if I've got the lead/lag function correct there (I've updated it based on a suggestion from @ram-rs), but it doesn't take too long to play round with all the combinations to get the right one.

FWIW, `%>%` (or `|>` now in native R) is a "pipe" operator that means "take the result of the last function and feed it into the first argument of the next function". It is used extensively in `tidyverse` scripts. More details here.
> 3 votes

---
Tags: r, genomics, bed, dataframe
---


# Question
Title: Nextflow: Properly chaining process outputs
I have written a pipeline that is composed of three processes. I seem to be misunderstanding how Nextflow passes the outputs from one process to the next or creates a queue of files. I am using a tumor and paired normal sample from SRA (SRX21185781) to prototype this pipeline.

My first process, `trim_QC` works fine and outputs the expected number of trimmed fastqs. However, the next process, `bwa_align` only creates a single bam from one of the two output from the preceding process. My terminal looks like this:

```
N E X T F L O W  ~  version 23.10.1
Launching `../main.nf` [sick_church] DSL2 - revision: beb03cc0c6

[6e/48d4f6] process > trim_QC (SRR25452858_normal)                  [100%] 2 of 2, cached: 2 ✔
[39/1b8cb2] process > bwa_align (SRR25452858_normal)                [100%] 1 of 1, cached: 1 ✔
[5f/e598c1] process > mapping_stats (SRR25452858_normal_sorted.bam) [100%] 1 of 1, cached: 1 ✔
...

```

and when running with `-with-trace` it is even more obvious only one bam file is being processed by the downstream processes:

```
task_id hash    native_id   name    status  exit    submit  duration    realtime    %cpu    peak_rss    peak_vmem   rchar   wchar
2   d3/0bd651   71127   trim_QC (SRR25452863_tumor) CACHED  0   2024-02-07 16:28:37.449 1h 21m 10s  1h 21m 10s  -   -   -   -   -
1   6e/48d4f6   71128   trim_QC (SRR25452858_normal)    CACHED  0   2024-02-07 16:28:37.459 1h 16m 32s  1h 16m 32s  -   -   -   -   -
3   39/1b8cb2   2754    bwa_align (SRR25452858_normal)  CACHED  0   2024-02-08 10:08:46.033 4h 18m 10s  4h 18m 10s  -   -   -   -   -
4   5f/e598c1   9282    mapping_stats (SRR25452858_normal_sorted.bam)   CACHED  0   2024-02-08 14:26:55.823 53.9s   53.8s   -   -   -   -   -

```

I am expecting for the `bwa-align` and `mapping_stats` processes to also be n of 2.

I have tried explicitly designating the `trim_QC.out` as a channel (`Channel.of(trim_QC.out)`) or using an operator to return a queue channel object (i.e., `.collect`, `.collectFile()`, etc.) but this either returns errors or fails.

If I try to visualize the outputs from the `trim_QC` process by adding this to my workflow: `trim_QC.out.trimmed_reads.view{it}`, it does show both file pairs in my terminal, but only the first one is processed by `bwa_align` and subsequent steps. Like so:

```
...
[6e/48d4f6] process > trim_QC (SRR25452858_normal)            [100%] 2 of 2, cached: 2 ✔
[39/1b8cb2] process > bwa_align (SRR25452858_normal)          [100%] 1 of 1, cached: 1 ✔
[5f/e598c1] process > mapping_stats (SRR25452858_normal_so... [100%] 1 of 1, cached: 1 ✔
hello

[SRR25452858_normal, /Users/me/Desktop/fastq2/work/6e/48d4f61cae83ff903abceaadf715e5/SRR25452858_normal_1_val_1.fq.gz, /Users/me/Desktop/fastq2/work/6e/48d4f61cae83ff903abceaadf715e5/SRR25452858_normal_2_val_2.fq.gz]
Analysis complete for SRR25452858_normal_1_val_1.fq.gz
Analysis complete for SRR25452858_normal_2_val_2.fq.gz

[SRR25452863_tumor, /Users/me/Desktop/fastq2/work/d3/0bd6515fcedfc060596b75dd90c28a/SRR25452863_tumor_1_val_1.fq.gz, /Users/me/Desktop/fastq2/work/d3/0bd6515fcedfc060596b75dd90c28a/SRR25452863_tumor_2_val_2.fq.gz]
Analysis complete for SRR25452863_tumor_1_val_1.fq.gz
Analysis complete for SRR25452863_tumor_2_val_2.fq.gz

```

I am unsure if I am misunderstanding what Nextflow is capable of or how channels work. What I am trying to accomplish is trim adapters and QC result, align trimmed reads followed by indexing, and finally calculate mappings statistics. Any advice is appreciated.

Here is my entire script:

```
#!/usr/bin/env nextflow
nextflow.enable.dsl = 2

params.fastq = "$projectDir/fastq2/*_{1,2}*"
params.index_dir = "$projectDir/RefGenome/BWAIndex"
params.ref = "hg38.fa"
params.out_dir = "$projectDir/Results"

workflow {

    fastq_ch = Channel.fromFilePairs(params.fastq)
    index_ch = Channel.fromPath(params.index_dir)
    ref_ch = Channel.of(params.ref)

    trim_QC(fastq_ch)
    bwa_align(index_ch, ref_ch, trim_QC.out.trimmed_reads)
    mapping_stats(bwa_align.out.sorted_bam)
}

process trim_QC {
    debug true
    tag "$sample_id"

    publishDir = [
                    [
                        path: { "${params.out_dir}/trimmed/" },
                        mode: 'copy',
                        overwrite: false
                    ],

                    [
                        path: { "${params.out_dir}/stats/" },
                        mode: 'copy',
                        pattern: "*fastqc*",
                        overwrite: false
                    ]
                ]

    input:
    tuple val(sample_id), path(fastq)

    output:
    tuple val("${sample_id}"), path("*_val_1.fq.gz"), path("*_val_2.fq.gz"), emit: trimmed_reads
    path("*fastqc*")

    script:
    def (read1, read2) = fastq
    """
    trim_galore --paired --fastqc ${read1} ${read2}
    """
}

process bwa_align {
    debug true
    tag "$sample_id"

    publishDir "${params.out_dir}/bwa_aligned/", mode: 'copy', overwrite: false

    input:
    path index_dir
    val ref
    tuple val(sample_id), path(read1), path(read2)

    output:
    path("${sample_id}_sorted.bam"), emit: sorted_bam
path("${sample_id}_sorted.bam.bai"), emit: bam_index

    script:
    """
    bwa mem -t 2 ${index_dir}/${ref} ${read1} ${read2} | samtools sort -@2 -o "${sample_id}_sorted.bam" -
samtools index ${sample_id}_sorted.bam > "${sample_id}_sorted.bam.bai"
    """
}

process mapping_stats {
    debug true
    tag "$sorted_bam"

    publishDir "${params.out_dir}/stats/", mode: 'copy', overwrite: false

    input:
    path(sorted_bam)

    output:
    path("${sorted_bam}_mapping_stats.txt")

    script:
    """
    samtools flagstat ${sorted_bam} > "${sorted_bam}_mapping_stats.txt"
    """
}

```

# Answer
> 2 votes
This is a result of nextflow channel types, queue- and value-channels. Nextflow will make fifo combinations of your channels to run `bwa_align` until either or both channel are empty. `index_ch` contains only one item, so once it is used there are no more valid combinations to start new `bwa_align` processes.

Either create your `index_ch` as value channel, or keep the index out of a channel entirely and refer to it as a parameter.

---
Tags: nextflow
---


title: Is ./1 a valid GT value to add to VCF?

text:
# Question
Title: Is ./1 a valid GT value to add to VCF?
I have a VCF file coming from a structural variant detection tool that doesn't add GT value. I can add GT values, for example by using this tool: https://vatools.readthedocs.io/en/latest/vcf\_genotype\_annotator.html. I know there exist the alternate allele (1), but I don't know if it's homozygous or heterozygous.

My collaborator wants a VCF with genotype. If an SV is called but not genotyped, I was wondering if ./1 could mean unknown allele and alternate allele for a diploid genome (working with human DNA). Although, if it is not a standard format, it may be a problem.

Is `./1` a valid value to add or does it have to be `0/1` or `1/1`? Thanks.

# Answer
> 1 votes
No, `./1` isn't valid. You can find the VCF specifications at https://samtools.github.io/hts-specs/. The current version is VCFv4.4 and this is the genotype section which hasn't really changed much from previous versions as far as I recall:

> GT (String): Genotype, encoded as allele value preceded by either of / or | depending on whether that allele is considered phased. The first phasing indicator may be omitted and is implicitly defined as / if any phasing indicators are / and | otherwise. The allele values are 0 for the reference allele (what is in the REF field), 1 for the first allele listed in ALT, 2 for the second allele list in ALT and so on. For diploid calls examples could be 0/1, 1 | 0, /0/1, or 1/2, etc. Haploid calls, e.g. on Y, male non-pseudoautosomal X, or mitochondria, should be indicated by having only one allele value. A triploid call might look like 0/0/1, and a partially phased triploid call could be |0/1/2 to indicate that the first allele is phased with another variant in the VCF. If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype). The meanings of the phasing indicators are as follows (see the PS and PSL fields below for more details on incorporating phasing information into the genotypes):
>
> * / : allele is unphased
> * | : allele is phased (according to the phase-set indicated in PS or PSL)
>
> For symbolic structural variant alleles, GT=0 indicates the absence of any of the ALT symbolic structural variants defined in the record. Implementer should note that merging a VCF record containing only symbolic structural variant ALT alleles with a record containing other alleles will result a change of the meaning of the GT=0 haplotypes from the record containing only symbolic SVs.

This is the bit you are most interested in:

> If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype).

So, `./1` doesn't mean anything. You could use `./.` which is the standard for "not genotyped", or you could simply not have a GT field at all since that isn't actually required by the standard, so a VCF file with no genotype field at all can be a perfectly valid and standard compliant VCF.

In other words, just think about what you really need. Do you absolutely need a GT field? If not, just leave it out and you're done. If you do, you can use `./.` if you want to say "I don't know what the genotype here is" (but be aware that some tools interpret `./.` as "this might not be a real call"), but don't use `./1` since that won't be understood by anyone.

---
Tags: vcf
---


# Question
Title: Bar Graph of Expression Data from Seurat Object
I am looking for a way to present some single-cell RNA sequencing data on expression of a certain gene. My dataset has 2 variables, cell type and condition. I'm looking to create a grouped aligned scatter plot bar graph of the expression of the gene, where the y-axis is expression levels, while X axis is the different cell types, with each cell type having 2 bars (1 for healthy, 1 for diseased). My dataset has 3 healthy and 3 diseased samples, but all of the data is integrated into a Seurat object. To first create an aligned scatter plot bar graph, what I did was generate a DotPlot for the expression of gene X in each sample, split by cell-type. Then I realized that DotPlot() prints out a data table with values including avg.exp, pct.exp and avg.scaled.exp. I used the avg.exp metric to plot the data, but I was posting this question to see if there are other ways I can access and present expression data from a Seurat object. Is the way I've presented accurate? I've posted a sample bar graph of what I'm looking for, and the main question is how I can gather the expression levels across cell types AND diseased or healthy condition.

Thanks!

# Answer
Regarding the other answer with the below code:

```
data.df <- object@meta.data[, c("celltypes", "condition")]
data.df["value"] <- object$RNA@data["geneX", ]
data.df <- data.df %>% mutate(group = paste(celltypes, condition, sep=":"))
data.df <- data.df %>% filter(value > 0)
data.df.2 <- data.df %>% group_by(group) %>% summarise(sd=sd(value),value=mean(value))
X <- str_split(data.df.2$'group', ":", 2, simplify=TRUE)

data.df.2$'celltypes' <- X[, 1]
data.df.2$'condition' <- X[, 2]

ggplot(data.df.2, aes(fill=condition, y=value, x=celltypes)) +
    geom_bar(color="black", position="dodge", stat="identity") +
    geom_errorbar(
        aes(ymin=value, ymax=value+sd), width=.2,
        position=position_dodge(.9)
    )

```

I find the plot matches more the results of Seurat differential expression if you DO NOT perform the step of removing the cells with 0 counts (`data.df <- data.df %>% filter(value > 0)`). To avoid having huge error bars you can use standard error of the mean instead of sd.
> 2 votes

# Answer
In single-cell data a bargraph with error bar is not usually shown since the gene expression does not have a normal distribution due to drop-outs. A violin-plot is commonly used. In your case, first set the active.ident to cell types

```
seuratObj <- SetIdent(seuratObj , value = "Cell.Types")

```

Then call the violin plot

```
VlnPlot(seuratObj , features = "geneX", split.by = "condition")

```

If you still specifically need bargraph, you can get the output of the violin plot as a dataframe and then draw the barplot.
> 1 votes

# Answer
Following the answer of @*Ajay*.

You can add `+ geom_boxplot()` after the command to turn the violin plot to a box plot.

As for the barplot, you can do it with these:

```
# the object is your seurat object
data.df <- object@meta.data[, c("celltypes", "condition")]
data.df["value"] <- object$RNA@data["geneX", ]
data.df <- data.df %>% mutate(group = paste(celltypes, condition, sep=":"))

# remove cells that do not express the gene
data.df <- data.df %>% filter(value > 0)

data.df.2 <- data.df %>% group_by(group) %>% summarise(sd=sd(value),value=mean(value))
X <- str_split(data.df.2$'group', ":", 2, simplify=TRUE)
data.df.2$'celltypes' <- X[, 1]
data.df.2$'condition' <- X[, 2]

ggplot(data.df.2, aes(fill=condition, y=value, x=celltypes)) +
    geom_bar(color="black", position="dodge", stat="identity") +
    geom_errorbar(
        aes(ymin=value, ymax=value+sd), width=.2,
        position=position_dodge(.9)
    )

```

Dropout is very usual in single-cell data. If you really want to show the expression with barplot, I would suggest removing the zero values or using denoised data for better visualisation.
> 1 votes

---
Tags: r, scrnaseq, seurat
---


title: Increasing memory usage in AWS instance

text:
# Question
Title: Increasing memory usage in AWS instance
I have 16 CPUs on my AWS EC2 instance and I am trying to use kSNP4 to generate an SNP matrix. I was able to successfully run this with over 150 samples - while it took a while to run it did complete - and now I want to run it with over 170 samples but when I use the same command it dies at the second step of the process. I have reason to believe it is due to the processors on my server that I am using more than the tool itself.

Does anyone know the best way to ensure maximum processing ability to run tools on an EC2 instance? Or what I can do to increase the power of the instance?

I know this is more than a technical question than a specific bioinformatic question but I would appreciate any insight!

# Answer
You need to diagnose for low memory and consider switching to an **R5** processor that prioritises memory. Basically it's easily solvable by switching hardware. There are certainly alot of other approaches but this is specialist HPC.

Personally, I would use Slurm ParallelCluster in AWS. This is very easy to set up from from the configuration file.

The overview of instance types here. https://aws.amazon.com/ec2/instance-types/
> 1 votes

---
Tags: linux, aws, cpu, ksnp4
---


title: Chromosome position -> ENSP

text:
# Question
Title: Chromosome position -> ENSP
I'd like a tabix file that, for every position in the genome, gives me a list of every ENSP within 5kb of that position.

This should be easy right... I know VEP /kinda/ gives me this, but I'd rather not run it for every postition in the genome

Can I use biomart to get the start/end position for every ENSP, then bedtools to somehow expand those coords?

# Answer
I used BioMart (and then tabix) to get the answer to this question:

```
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

BiocManager::install("biomaRt")

## F**kery:
install.packages("devtools")
devtools::install_version("dbplyr", version = "2.3.4")

library(biomaRt)

listEnsembl()

ensembl <- useEnsembl(
  biomart = "ensembl",
  dataset = "hsapiens_gene_ensembl",
  GRCh = 37
)

chrs <- c(1:22, "X", "Y", "MT")

genes <-
  getBM(
    attributes = c(
      "chromosome_name",
      "start_position",
      "end_position",
      "ensembl_gene_id",
      "ensembl_transcript_id",
      "ensembl_peptide_id",
      "hgnc_symbol"
    ),
    filters = "chromosome_name",
    values = chrs,
    mart = ensembl
  )

```
> 1 votes

---
Tags: ensembl, bed, bedtools, biomart
---


# Question
Title: Salmon Multiple File SCRIPT giving error
Hi I am using Salmon quantitation for multiple fastq paired files using following code

```
salmon_index="human_salmon_index"
fastq_dir="fastq/"
for dir in "${fastq_dir}"/SRR*
do
r1_file=$(find "$dir" -name "*_1.fastq")
r2_file=$(find "$dir" -name "*_2.fastq")
samp=$(basename "$dir")
echo "Processing sample ${samp}"
salmon quant -i "$salmon_index" -l A -1 "$r1_file" -2 "$r2_file" -p 10 -o "salmon_out/${samp}_quant"
done

```

and getting this error after it iterates over all files and start initial processing steps.

```
Processing sample SRR27756943_2.fastq
Version Info: This is the most recent version of salmon.
### salmon (selective-alignment-based) v1.10.0

Exception : [
The following errors were detected with the read files

======================================================

ERROR: file [] does not appear to exist!

```

Any solution?

# Answer
> 0 votes
Are you sure the `find` operations return exactly one file each? The error seems to state that there were no files found for some run. Or it could also be that you have a file that matched your glob pattern that's not a dir, so the `find` operation fails.

Try this:

```
salmon_index="human_salmon_index"
fastq_dir=fastq
for dir in ${fastq_dir}/SRR*/
do
r1_file=$(find $dir -name "*_1.fastq")
r2_file=$(find $dir -name "*_2.fastq")
samp=$(basename $dir)
echo "Processing sample ${samp}"
echo salmon quant -i "$salmon_index" -l A -1 "$r1_file" -2 "$r2_file" -p 10 -o "salmon_out/${samp}_quant" | tee salmon_commands.cmd
done

```

Verify the contents of `salmon_commands.cmd` then use `bash` or `parallel` to run the actual commands.

---
Tags: rna-seq, sequence-alignment, bash, salmon, quantification
---


# Question
Title: After running nf-core, is there a way to map a gene_id to a specific gene's DNA sequence?
I have been running nf-core in Python and it works great! But I have a seemingly simple question that I'm struggling to find an answer for online. After running the nf-core pipeline on my RNA sequencing data, Salmon is able to provide gene counts in a file named `salmon.merged.gene_counts.tsv`. This file has essentially all the information I am seeking with all counts associated to a unique gene\_id. But, for a large number of these gene\_ids where no specific gene\_name is known, I would like to further find the corresponding DNA sequence (or corresponding transcript's RNA sequence) for each gene\_id. Accordingly, my question is: where can I find the information to map the gene\_id to the corresponding nucleotides of the gene's (or original transcript's) sequence?

(Apologies in advance if my question is missing information or overly simply and I'm overlooking something obvious — I am quite new to this all.)

# Answer
#### fishing out tp53 RNA-Seq reads

* go to NCBI Taxonomy Browser
* put `Anas platyrhynchos` as a search term
* click on `Anas platyrhynchos`
* on the right you will have a table with `SRA Experiments`, click on `11,679`
* restrict to `RNA` and `paired` on the left
* send it to `Run Selector`
* order it by the `Bases`
* select few top runs, like `SRR16082819`

Now it depends what is your exact goal. You may (separate approaches):

* get mallard transcriptome as complete as possible, use `kallisto` just to check that a given pair of FASTQs contains bunch of reads derived from tp53
* use `STAR` to align it to the genome with the contig containing TP53 gene
* one can also try to fish out the reads out of FASTQ using tp53 cDNA sequence, but this will likely miss i.e. putative alternative 5- & 3-prime part of the gene
> 0 votes

---
Tags: rna-seq, gene-expression, transcriptome, star, salmon
---