cannin/biostars_qa huggingface dataset:
Question: I'm a Computer Science undergraduate writing my thesis which is about analyzing open source bioinformatics projects and extracting standalone modules. I need those projects to be written in Java and preferably be stable or mature. So far I've found a bunch of projects via SourceForge and I was wondering if there are any notable ones that are not hosted there. A different repository that has a bioinformatics category would be helpful too.

Here's a list of the most promising projects I have so far:

- Jenetics
- JGap
- Jmol
- Juicebox
- Picard
- RDP Classifier
- Shim
- Tassel
- The Chemistry Development Kit
- Toxfree
- TreeView
- VarScan

Thanks in advance.

Answer: - [BBMap][1]
- [FastQC][2]

[1]: https://sourceforge.net/projects/bbmap/
[2]: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Question: Is there a way to output a gene expression matrix corrected for batch (I know I can use ComBat) and surrogate variables from SVA? I would like to visualize the results of a limma differential expression however they are confounded and thus not amenable to visualization.

Answer: If you have:

- `y` as the gene expresion matrix
- `mod` as the model matrix you sent to sva (the full model)
- `svs` as svobj$sv where svobj is the output from the sva function

then you can use the function below to get a "cleaned" version of the matrix with the surrogate variables removed.

```
cleanY = function(y, mod, svs) {
X = cbind(mod, svs)
Hat = solve(t(X) %*% X) %*% t(X)
beta = (Hat %*% t(y))
rm(Hat)
gc()
P = ncol(mod)
return(y - t(as.matrix(X[,-c(1:P)]) %*% beta[-c(1:P),]))
}
```

I modified this from a post by Andrew Jaffe [here][1].

[1]: http://permalink.gmane.org/gmane.science.biology.informatics.conductor/42857
Question: Hi guys,

So I've been trying to `bgzip` around 100 VCF files in parallel, but although the jobs are submitted and files get created there's definitely something wrong.

So far I've been trying:

$ for file in *.vcf
> do
> bsub /foo/bar/bgzip $file
>> $file.gz

What is the correct way to do this?

Thanks in advance!

Answer: Use [GNU parallel][1]

parallel bgzip {} ::: *.vcf

or `xargs`:

ls *.vcf | xargs -P 10 bgzip

[1]: http://www.gnu.org/software/parallel/
Question: Hi all,

I am interested in calling SNPs for some alignment data that I have. I am looking at solutions with `samtools mpileup` and `freebayes`.

I noticed that freebayes has an option to specify how many observations have to support a variant before it is called (the `-C` option).

Is there something similar for the `samtools mpileup` option?

I've been using `samtools` as follows:

`samtools mpileup -uf reference.fasta Sample_sorted.bam -d 8000 | bcftools call -mv -Oz -o variants.vcf`

However, I don't know how many observations supported the resulting variants.

An example SNP call from the above command is as follows:

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	Sample_sorted.bam
Ecoli_genome	5159	.	A	G	228	.	DP=303;VDB=0.658475;SGB=-0.693147;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;AC=2;AN=2;DP4=0,1,120,171;MQ=60	GT:PL	1/1:255,255,0
Ecoli_genome	5162	.	A	T	228	.	DP=310;VDB=0.633206;SGB=-0.693147;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,120,172;MQ=60	GT:PL	1/1:255,255,0
...

Any insight is appreciated!

Answer: "DP4=0,1,120,171" tells you the number of reads at that position, in this format:

DP4=refPlus,refMinus,varPlus,varMinus

So, in that case (for the first variant) there was 1 read supporting the reference, and 120+171=291 supporting the variant.

Incidentally, if you're trying out variant callers, the BBMap package also has one (callvariants.sh) that is extremely fast and easy to use. It has a flag "minreads" that can be set to specify the number of observations required to call a variant, as well as various other filters.
Question: Hi all,
I've been given a set of **RNA-Seq** fastqs (4 samples) and a collaborator would like to know how many defined **intron-exon** junction in his favorite gene are detected in those reads.

My idea was to clean the reads with cutadapt, align with Bowtie and create a custom program to loop over the cigar string and detecting the junctions.

or is there an existing tool for this ?
Bonus if the tool can tell me wether two junctions belong to the same transcript.

Answer: If you use the STAR aligner, it will output `SJ.out.tab` with high confidence collapsed splice junctions in tab-delimited format.
Question: What is the best tool for extracting copy number alterations from a .maf file? I would like to generate a .cna file which only contains inserts and deletions.

Thanks.

Answer: You cannot do this due to the fact that the MAF [format] does not encode any copy number information specifically. Insertions and Deletions (indels) are regarded in this context as mutations, not copy number alterations. If you literally want to just extract indels from your MAF data, then that should be easy by simple `subset()` commands (in R).

There is a brief tutorial here about how one can integrate MAF data with copy number data: https://www.bioconductor.org/packages/devel/bioc/vignettes/maftools/inst/doc/maftools.html

If something else, then please elaborate further on what you are aiming to do.

Kevin
Question: I am doing a DESeq2 comparison with different levels and one factor. To do this, I have performed the analysis in two different ways.

First, putting all the samples in the same DESeq object and then extracting each comparison:

> sampleinfo
FileName	SampleName	Status
A_1_count	A_1	A
A_2_count	A_2	A
B_3_count	B_3	B
B_4_count	B_4	B
C_5_count	C_5	C
C_6_count	C_6	C
D_7_count	D_7	D
D_8_count	D_8	D
E_9_count	E_9	E
E_10_count	E_10	E

dds <- DESeqDataSetFromMatrix(countData = cts,
colData = sampleinfo,
design = ~ Status)

dds$Status <- relevel(dds$Status, ref = "E")

And the results:

dds <- DESeq(dds)
res_A <- results(dds,name="Status_A_vs_E")
res_B <- results(dds,name="Status_B_vs_E")
res_C <- results(dds,name="Status_C_vs_E")
res_D <- results(dds,name="Status_D_vs_E")

And doing these comparisons one by one separately on different DESeq objects.

> sampleinfo_A
FileName	SampleName	Status
A_1_count	A_1	A
A_2_count	A_2	A
E_9_count	E_9	E
E_10_count	E_10	E

> sampleinfo_B
FileName	SampleName	Status
B_3_count	B_3	B
B_4_count	B_4	B
E_9_count	E_9	E
E_10_count	E_10	E

dds_A <- DESeqDataSetFromMatrix(countData = cts_A,
colData = sampleinfo_A,
design = ~ Status)

dds_B <- DESeqDataSetFromMatrix(countData = cts_B,
colData = sampleinfo_B,
design = ~ Status)

And the results:

dds_A <- DESeq(dds_A)
res_A <- results(dds_A)

dds_B <- DESeq(dds_B)
res_B <- results(dds_B)

> (Repeat for each condition)

However, the results give me different between the 2 methods. Does anyone know why is this happening? How it is the correct way to compare all to E?

Thank you!

Answer: > However, the results give me different between the 2 methods. Does
> anyone know why is this happening?

Because of the dispersion estimates. By including more samples you will get better dispersion estimates
Question: Hi,

I'm trying do some permutation testing on microarray data in R and the sample labels look like this:

labels <- c("case_1", "case_2", "case_3", "case_4", "case_5",
"control_1", "control_2", "control_3", "control_4", "control_5")

To permute them I realize I could just repeatedly use "sample(labels)" to do random permutations, but I'd like to look ALL possible permutations.

I think there should be ((10 choose 5) / 2) = 126 possibilities.

Does anyone know if there's an easy way to generate these though?

Thanks!

Answer: Using the standard R libaries, you can do it with the combn() function

> a = c(1,2,3,4)
> combn(a,3)
[,1] [,2] [,3] [,4]
[1,] 1 1 1 2
[2,] 2 2 3 3
[3,] 3 4 4 4

Your combinations are in the columns.
Question: Hi!

I have a pile of unindexed reads, including some PhiX, from a poor run on a NexSeq. Can anyone suggest some tools or techniques I can apply to this data just to try get a better idea of why the run went poorly?

Much appreciated.

Answer: Besides the obvious issue of contamination (which you can check by @WouterDeCoster's method) you should also check on the following.

NOTE: If the reads are truly "un-indexed" (as in having "no index read" as opposed to having landed in the "undetermined" file since they can't be de-multiplexed using indexes expected to be present in the run, then the following does **not** apply).

If you are using bcl2fastq v.2.17.* and have access to the flowcell folder you could look at the `DemuxSummaryF1L*.txt` files in folder containing demux data. These files list all "undetermined" barcodes present (in addition to the ones properly demultiplexed) in that specific lane (there are other ways of counting the barcodes should you only have access to the undetermined*.fastq.gz file).

One common issue is poor selection of indexes for pooling. Having low diversity indexes along with a borderline overloaded run can lead to "N" basecalls in indexes leading to problems with demultiplexing. You can overcome this partially by allowing one or two mismatches in tag reads but depending on the selection of indexes/length of your tag read this may or may not work well.

Illumina has published guides on what indexes to use for pooling for variable smaller number of samples (for both nextera/truseq kits) if that indeed is the issue.
Question: I asked this question a few days ago but did not get any answer. I tried finding it on my but could not really get it. I have orthofinder results 'Gene_count' which gives gene count in each orthogroup in each genome. Its format is like;
>OG0000032	10	11	1	13	0	35

I also have SpeciesTree_rooted that is generated by orthofinder. I need to know if someone here has used orthofinder results and some modifications (if they did) to get successful gene families expansion/contraction results from CAFE or any other software please? I am really looking for an answer!

Answer: I have been looking into something similar (maybe, not sure myself) and I didn't know about cafe. It looks good, thanks for bringing it up. The following R code should adapt orthofinder 2.x output to [cafe5](https://github.com/hahnlab/CAFE5).

Lots of disclaimers: This code works for me in that it makes cafe complete without errors. Whether it makes something sensible I'm not sure as I'm still getting to grips with the whole idea of expansion/contraction of gene families.

First, make orthofinder's species tree ultrametric. It should be already binary and rooted as required by cafe.

```
library(ape)
library(data.table)

tre <- read.tree('Species_Tree/SpeciesTree_rooted.txt')
stopifnot(is.binary(tre))
stopifnot(is.rooted(tre))

if(is.ultrametric(tre)) {
utre <- tre
} else{
utre <- chronos(tre)
}
write.tree(utre, 'cafe/SpeciesTree_rooted_ultra.txt')
```

Then prepare the table of counts. I use the N0 output since the original orthogroups are deprecated in favour of the phylogenetic hierarchical orthogroups:

```
hog <- fread('Phylogenetic_Hierarchical_Orthogroups/N0.tsv')
hog[, OG := NULL]
hog[, `Gene Tree Parent Clade` := NULL]
hog <- melt(hog, id.vars='HOG', variable.name='species', value.name='pid')
hog <- hog[pid != '']
hog$n <- sapply(hog$pid, function(x) length(strsplit(x, ', ')[[1]]))

# Exclude HOGs with lots of genes in a one or more species.
# See also cafe tutorial about filtering gene families
keep <- hog[, list(n_max=max(n)), HOG][n_max < 100]$HOG
hog <- hog[HOG %in% keep]

# Exclude HOGs present in only 1 species
keep <- hog[, .N, HOG][N > 1]$HOG
hog <- hog[HOG %in% keep]

counts <- dcast(hog, HOG ~ species, value.var='n', fill=0)
counts[, Desc := 'n/a']
setcolorder(counts, 'Desc')
fwrite(counts, 'hog_gene_counts.tsv', sep='\t')
```

Run cafe5 with your favorite settings:

```
nohup cafe5 -i cafe/hog_gene_counts.tsv -t cafe/SpeciesTree_rooted_ultra.txt -o cafe/results --cores 30 &
nohup cafe5 -i cafe/hog_gene_counts.tsv -t cafe/SpeciesTree_rooted_ultra.txt -o cafe/results_k3 --cores 30 -k 3 &
nohup cafe5 -i cafe/hog_gene_counts.tsv -t cafe/SpeciesTree_rooted_ultra.txt -o cafe/results_p_k3 --cores 30 -p -k 3 &
```
Question: I got a vcf file (header and individual ID have been removed form the vcf file), which contains only genotype information for 1000 positions, looks like this (for the fist position)

head all_unique.vcf

5	49444858	rs12513792	C	A	100	PASS	AC=3353;AF=0.669529;AN=5008;NS=2504;DP=17854;EAS_AF=0.5625;AMR_AF=0.7075;AFR_AF=0.7489;EUR_AF=0.672;SAS_AF=0.6431;AA=.|||;VT=SNP	GT	0|1	1|0	0|1	1|0	0|1	0|1	1|1	1|0	1|0	1|1	1|1	1|0	0|1	1|1	1|1	0|1	1|1	0|1	0|1	0|1	0|0	0|0	0|1	0|1	1|1	0|1	1|0	1|1	1|0	1|0	1|0	0|0	1|1	1|0	1|1	1|0	1|1	1|1	1|1	1|1	1|1	0|0	1|1	0|1	0|1	1|1	1|1	1|0	1|1	1|1	0|1	1|1	1|1	1|1	1|1	1|0	1|1	0|1	0|0	1|1	1|0	1|1	1|1	1|1	1|1	1|1	1|0	1|0	1|1	0|1	0|1	1|1	1|0	1|0	1|0	0|0	0|1	1|0	0|1	1|1	0|0	1|0	1|0	0|1	1|1	1|1	1|1	1|1	1|0	1|1	1|1	1|0	0|1	1|1	0|1	0|1	0|1	0|0	1|1	0|1	1|1	0|1	0|1	1|1	0|1

I have another file which contains the individual ID:

head individual_id

HG00096
HG00097
HG00099
HG00100
HG00101
HG00102
HG00103
HG00105
HG00106
HG00107
HG00108
HG00109

the first individual in the ID file corresponds to the first genotype in the vcf file, and so on. For each position,I want to attach ID to each genotype and get an output and save the output in a seperate file like this :

head desired_output_pos49444858
HG00096 0|1
HG00097 1|0
HG00099 0|1
HG00100 1|0
HG00101 0|1

Does anyone have an idea how can I make this file?
Thanks

Answer: assuming that the number of lines in individual_id is the same as the number of genotypes in the vcf...

cat all_unique.vcf | while read L ; do P=$(echo "$L" | cut -f2);paste <(cat individual_id ) <(echo "$L"| cut -f 10- | tr "\t" "\n" ) > "desired_output_pos${P}";done
Question: I ran HISAT2 (index built using a transcriptome multi fasta) intending that it won't perform gapped alignment.
I use following script to run HISAT:

INDEX=./indices/hisat/transcriptome
FASTQ=$1
OUTPUT=./transcriptome_aligned/$2.sam
./software/hisat-0.1.6-beta/hisat \
-q \
-p 2 \
--no-spliced-alignment \
--end-to-end \
-x $INDEX \
-U $FASTQ \
-S $OUTPUT

Should I still expect gapped alignment in my SAM file?
I have records like this in the SAM output.

SRR2144041.255	0	YCL025C	274	255	16M1I33M	*	0	0	CAGGCTCAAGAACTAGAAAAAAAATGAAAGTTCGGACAACATAGGCGCTA	CCCFFFFFHHHHHJJJJJJJJJJJJJJJIJGIIIIJJJJJIJJIIJJJHH	AS:i:-8	XN:i:0	XM:i:0	XO:i:1	XG:i:1	NM:i:1	MD:Z:49	YT:Z:UU	NH:i:1

This shows that HISAT2 is still performing gapped alignment even with `--end-to-end` and `--no-splice-alignment` parameters.

I'm trying to use the output SAM for `rsem-calculate-expression` but it returns following error due to presence of gapped alignment:

rsem-parse-alignments ./indices/rsem/rsem ./rsem_output/sample.temp/sample ./rsem_output/sample.stat/sample ./transcriptome_aligned/sample.bam 1 -tag XM
Read SRR2144041.836747: RSEM currently does not support gapped alignments, sorry!

"rsem-parse-alignments ./indices/rsem/rsem ./rsem_output/sample.temp/sample ./rsem_output/sample.stat/sample ./transcriptome_aligned/sample.bam 1 -tag XM" failed! Plase check if you provide correct parameters/options for the pipeline!

How do I make sure that HISAT2 doesn't perform gapped alignment? Should I filter the output for using `grep -v XO:i:0`?

**EDIT:**

I checked RSEM manual and found that in order to avoid gapped alignments using Bowtie2, RSEM uses following Bowtie2 parameters:

--sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1

I wonder what is the equivalent of `--gbar` in HISAT2

Thanks

Answer: Solved it:

Included following two parameters in HISAT2 run command to make gap opening and extension insanely difficult.

--rdg 10000,10000
--rfg 10000,10000

My HISAT2 script now looks like this:

INDEX=./indices/hisat/transcriptome
FASTQ=$1
OUTPUT=./transcriptome_aligned/$2.sam
./software/hisat-0.1.6-beta/hisat \
-q \
-p 2 \
--no-spliced-alignment \
--end-to-end \
--rdg 10000,10000 \
--rfg 10000,10000 \
-x $INDEX \
-U $FASTQ \
-S $OUTPUT

Hope it helps, if anybody is trying to run RSEM downstream of HISAT2.

This solution might not be foolproof. Improvements on this answer are invited.
Question: Hi,

I'm trying to reformat my fasta file headers in terminal on my mac. Currently they look like this:

```
tpg|Magnaporthiopsis_incrustans|JF414846
tpg|Pyricularia_pennisetigena|AB818016
tpg|Inocybe_sororia|EU525947
```

However, I need them to look like this, with the species last:

```
tpg|JF414846|Magnaporthiopsis_incrustans
tpg|AB818016|Pyricularia_pennisetigena
tpg|EU525947|Inocybe_sororia
```

Could anyone advise on how best to do this?

Thanks in advance!

Answer: Heres a slightly dodgy bash only script, if you don't want to/can't install BioAwk/Perl/Python etc. It assumes the input files are perfectly formatted fasta's with headers as per the OP, i.e. `>xxxx|yyyyy|zzzzzzzzzz`

(I haven't tested it with fasta's that weren't already linearised but it should work...)

#!/bin/bash
# $ script.sh sequences.fasta

while read line ; do
if [ ${line:0:1} == ">" ] ; then
desc1=$(echo $line | cut -d '|' -f 1)
desc2=$(echo $line | cut -d '|' -f 2)
desc3=$(echo $line | cut -d '|' -f 3)
else
seq=$line
echo -e "$desc1|$desc3|$desc2"
echo -e "$seq"
fi
done < $1

In this script, the sequence fasta is looped over line by line, if it hits a header line (i.e. starts with a ">"), which in `bash` is expressed here as `${line:0:1}` (the zero-th to first character range, AKA, the first character, is a ">"), the string is split on the `|` (pipe symbol) delimiter, and `cut` then picks out each section of the header and assigns it to a variable.

If the string didn't match the ">", it must be a sequence line, and can just be re-emitted as it is.

`echo` then prints each variable, separated by `|` again, in the order you want (1, 2, 3 becomes 1, 3, 2).

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

EDIT: a slightly neater version

#!/bin/bash
# $ script.sh sequences.fasta

while read line ; do
if [ ${line:0:1} == ">" ] ; then
IFS='|' read -a header <<< "$line"
else
seq="$line"
echo -e "${header[0]}"\|"${header[2]}"\|"${header[1]}""\n""$seq"
fi
done < $1

Same process basically, except this time we can make use of the inbuilt "IFS" (internal field separator), and tell it to use `|` instead. In conjuction with `read -a` we can turn a string, in this case the header line, in to an array, saving us a few lines of code.

We also don't have to handle multiple variable names explicity, since bash array syntax lets you address individual items with `"${array[1]}"` for example (which is zero based, so is actually the second element). `"${array[@]}"` would refer to every element.

If we include a `"\n"`, we can also emit the fasta in a single echo line, though I personally would argue this is less neat (certainly less readable).
Question: Hi guys,

As far as I've read, both `SomaticSniper` and `Varscan` variant calling algorithms have been developed by the same group. Do you have any experience of which one is more accurate? Or in which conditions should I use one or the other?

Thank you in advance,

Answer: It just depends on what you're doing. Somatic Sniper has largely been superceded by a new generation of callers (Mutect, Strelka, etc) but is still useful in the niche case of having tumor contamination in your normal sample (very few other callers are designed to handle this case). In general, VarScan has kept up with the times pretty well.

That said, relying on a single caller is probably not your best bet. Each has strengths and weaknesses, and an ensemble approach followed by intelligent merging is probably your best bet for maximizing your performance.
Question: Hi,
I am looking for a better explanation for N50 in genome assembly. As per my understanding, N50 is the length of the contigs which covers 50% of genome. Am I right ?
Also, say for example, if I have 2 tools which give N50 as 500 and 1000 respectively, which of these would be a better tool ?
Thanks.

Answer: Your definition/understanding of N50 is somewhat correct indeed.

the way you calculate N50 is : you order your contigs large to small, then you start making the cumulative sum of the lengths of the contigs until you have >50% of your assembly , that number of contigs you need to get to 50% is the L50, and it's length is N50.

Intuitively one should go for the assembly with the highest N50 (1000 in this case), but N50 alone is not a good measure of performance, also total assembled size etc are of importance (NG50 might help here a little).
Question: I am using the -num_threads option for magicBLAST and running from a linux shell. However only one thread is used (tried on 2 machines). The file.fa is only 20MB.

magicblast -query file.fa -db nt -num_threads 16 -out out.sam -no_unaligned

Any ideas?

Answer: After 20 min was using multiple threads. I presume there is a stage when magicBLAST only uses a single thread when starting. NB I also reduced num_theads to 4 but havent tested if that has anything to do with it (24 thread machine)
Question: Is there gold standard data for ancestral reconstruction? I am asking because I am trying to create a new approach for ancestral (gene/protein sequence) reconstruction (using a algorithm that has not been used before for this) and I wanted to compare the performance of current tools like MrBayes, MEGA, Mesquite etc.

Is there a way to compare their performances based on some gold standard data? I would like to use the same data to check the efficiency of my new approach.

It's hard for me to believe that there wouldn't be any such dataset, given that so many tools provide this service.

Any help would be appreciated.

Answer: For this, I would use simulation, see for example [Computational reconstruction of ancestral DNA sequences.][1]
If you want experimental data, have a look at [An experimental phylogeny to benchmark ancestral sequence reconstruction][2].

[1]: https://www.ncbi.nlm.nih.gov/pubmed/18629667
[2]: https://www.nature.com/articles/ncomms12847
Question: Dear all,

I'm new in fungal ITS (internal transcribed spacer) metabarcoding analysis, so please ask if you don't understand my question.

Our lab received Illumina MiSeq sequencing result of our 150 samples from a company, but they give us only the R1 and R2 fastq files, without the consensus sequences. Our lab used to have 454 sequencing previously, and usually we get the consensus sequences. My question is, should we assemble those R1 and R2 into a consensus sequence before further analysis?

We have a pipeline for further analysis which is based on the 454 data, and I'm afraid that if we don't have consensus sequences beforehand, the R1 and R2 sequences would be recognized as two different ITS by the pipeline.

Any idea?

And thanks a lot in advance.

Answer: After a few months, I've learnt one or two things about ITS analysis, and one of them is that R1 and R2 reads might be used without merging them beforehand. In fact, we might fail to identify many species using only the merged paired-reads (it's pointed out in [this paper][1], for example).

[1]: http://dx.doi.org/10.1111/nph.12923
Question: Hello,

I want to generate a gene arrow map showing all the details like gene features, annotation etc. I have the data in the following format:

start	stop	strand	function
208	885	+	Chromosomal replication initiator protein DnaA
1576	1157	-	FIG00958722: hypothetical protein
2174	1578	-	Flavoprotein WrbA
2524	2171	-	Arsenate reductase (EC 1.20.4.1)
2649	3884	+	Inner membrane protein YihY, formerly thought to be RNase BN
4034	4228	+	FIG00956556: hypothetical protein

Please let me know about any software that can convert the above mentioned information to well formatted gene arrow map.

Thanks

Answer: using awk :-D

awk -F '\t ' -v chromStart=208.0 -v chromEnd=4228.0 '/^start/{next;} {L=60.0;x1=L*(int($1)-chromStart)/(chromEnd-chromStart);x2=L*(int($2)-chromStart)/(chromEnd-chromStart);if(x1>x2) { t=x1;x1=x2;x2=t;} for(i=0;i=int(x1) && i<=int(x2)) { c=$3=="+"?">":"<" } ; printf("%c",c);}printf(" %s\n",$4);}' input.tsv

>>>>>>>>>>> Chromosomal replication initiator protein DnaA
<<<<<<< FIG00958722: hypothetical protein
<<<<<<<<<< Flavoprotein WrbA
<<<<<< Arsenate reductase (EC 1.20.4.1)
>>>>>>>>>>>>>>>>>>> Inner membrane protein YihY, formerly thought to be RNase BN
>>> FIG00956556: hypothetical protein

else see

* https://www.biostars.org/p/2398/
* https://www.biostars.org/p/80956/
* https://www.biostars.org/p/93448/
* ....