cannin/biostars_qa huggingface dataset:
Question: Hello,

I would like to understand what is the best way of using bwa in parallel in a SLURM cluster. Obviously, this will depend on the computational limits that I have as user.

bwa software has an argument "-t" specifying the number of threads. Let's imagine that I use `bwa mem -t 3 ref.fa sampleA.fq.gz`, this will mean that bwa split the job on three tasks/threads. In other words, it will align three reads at a time in parallel (I guess).

Now, if I want to run this command on several samples and in a SLURM cluster, Shall I specify the number of tasks as for bwa mem, and specify the number of CPUs per task(for instance 2)? Which would be:

sbatch -c 2 -n 3 bwa.sh

where bwa.sh containes:

cat data.info | while read indv; do
bwa mem -t 3 ref.fa sample${indv}.fq.gz
done

Do you have any suggestion? Or can you improve/correct my reasoning?

Answer: One good idea is to use a job array that takes a samplesheet in input. Each line specifies one sample. Then using SLURM_ARRAY_TASK_ID variable in the slurm bash script you can acces the respective line in your sample sheet and execute a bwa mem instance with the number of CPU specified in the slurm parameters.

something like for 5 CPU - 20G RAM and 10 samples. You may want to tune the size of the RAM required depending on the reference genome.

#!/bin/bash
#SBATCH --nodes=1 --ntasks-per-node=5 --mem-per-cpu=20G
#SBATCH --array=1-10
#SBATCH -J jobName

samplesheet="sample_sheet.txt"

threads=$SLURM_JOB_CPUS_PER_NODE

r1=`sed -n "$SLURM_ARRAY_TASK_ID"p $samplesheet | awk '{print $1}'` # fastq.gz file

bwa mem -t $threads ref.fa $r1

very simple sample_sheet.txt :

sample_1_R1.fastq.gz
sample_2_R1.fastq.gz
sample_3_R1.fastq.gz
sample_4_R1.fastq.gz

You may change it something more detailed. Here for example it's paired-end data. And column 1 contains the name of the sample.

sample_1	sample_1_R1.fastq.gz	sample_1_R2.fastq.gz
sample_2	sample_2_R1.fastq.gz	sample_2_R2.fastq.gz
sample_3	sample_3_R1.fastq.gz	sample_3_R2.fastq.gz
sample_4	sample_4_R1.fastq.gz	sample_4_R2.fastq.gz

Then within the script you can change the section :

r1=`sed -n "$SLURM_ARRAY_TASK_ID"p $samplesheet | awk '{print $1}'` # fastq.gz file

to

name=`sed -n "$SLURM_ARRAY_TASK_ID"p $samplesheet | awk '{print $1}'`
r1=`sed -n "$SLURM_ARRAY_TASK_ID"p $samplesheet | awk '{print $2}'`
r2=`sed -n "$SLURM_ARRAY_TASK_ID"p $samplesheet | awk '{print $3}'`
Question: Hi -

So I'm just wondering how exactly the BWT works in creating "runs" of the same character in the last column after sorting the cycled string. Is it because of some frequency of nucleotides occurring conditionally that I'm not aware of? I get that in the English language I'm pretty sure that there are cases where a character is more likely to appear after another character (that's how BWT's effectiveness is explained in all texts i've found so far), but does this also apply to nucleotides?

Answer: You are right, the BWT will have a hard time compressing random sequences of nucleotides as (uniform) random data is by definition hard to compress. However, genetic sequences are far from random. Just think of GC-bias, codon-bias, motifs, patterns in various forms (promotors, TATA-box), duplications, …. All of these reduce the "randomness" (entropy) of the data and instead increase the repetitiveness which the BWT then can exploit.
Question: Given a de novo assembly and a reference assembly, what methods have you tried / would you recommend for determining structural variations?

Answer: For this paper [An integrated map of structural variation in 2,504 human genomes][1] my tiny part was to validate complex structural variation in long read TruSeq data.

It was about 3 years ago so I'm not sure if there are better approaches. But we created breakpoint contigs across putative SV breakpoints using Velvet.

Then I took the breakpoint contigs (there are many possible ones generated by Velvet) and I used BLAT to align them to the reference genome.

Using the BLAT results I was able to parse out the precise breakpoints.

Like I said it's rather labor intensive and I'm sure there's a better way of doing it. But this might be a good lead!

[1]: http://www.nature.com/nature/journal/v526/n7571/abs/nature15394.html
Question: I have been using the DESeq VST method on gene counts produced by Htseq-count as follows:

cds <- newCountDataSet(countData = dat, conditions = factor(conditions))
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds, sharingMode = "gene-est-only", method = "pooled", fitType = "local")
vst <- getVarianceStabilizedData(cds)

But honestly, I do not understand what exactly the getVarianceStabilizedData() function does. Can someone explain in simple terms:

1. Why is it necessary to normalize raw count data? Why can't we use the raw count data?
2. How exactly are we normalizing the raw count data using `getVarianceStabilizedData()` function?
3. Should the conditions parameter in the newCountDataSet() function match the conditions between which you want to find differentially expressed genes? For e.g. I have both cases & controls as well as males & females. So should I include both the information in the conditions parameter or just cases & controls?

I know these questions can be searched for easily, and I did. But I want a simple explanation from someone who uses these methods regularly to clear my concepts.

Answer: The vignettes in DESeq2 (which you should prefer using these days) describe these things and why you'd want to use them, under these sections:

- http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.pdf Under "Data transformations and visualization; and
- http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/beginner.pdf Under "working with rlog-transformed data"

The main image you want to have in your mind is the one in Figure 3 of the first link. The point of these transforms is to reduce (ideally eliminate) dependence of the variance on the mean. The second link above has this paragraph which sums things up quite nicely:

> Many common statistical methods for exploratory analysis of multidimensional data, especially methods for clustering and ordination (e.g., principal-component analysis and the like), work best for (at least approximately) homoskedastic data; this means that the variance of an observable quantity (i.e., here, the expression strength of a gene) does not depend on the mean. In RNA-Seq data, however, variance grows with the mean.

You should prefer to use these when you are doing downstream analysis on your count data that doesn't involve testing for differential expression using the statistical methods developed for count data. These scenarios include doing things like clustering, or PCA over your expression data or using the data as input to another machine learning algorithm.

Take the time to read the two vignettes above, as well as the [DESeq2 preprint][1] to get a better understanding of this (and many other things related to differential expression analysis with this software) ... the authors have gone to great lengths to document their software and methodology quite thoroughly.

[1]: http://biorxiv.org/content/biorxiv/early/2014/02/19/002832.full.pdf
Question: I'm using TopGo to do GO analyses. When I get my enriched terms, the long term names get corrupted. ```I have some text then ...```. For example:

```
GO:0000122	negative regulation of transcription by ...
GO:0006886	intracellular protein transport
GO:0010976	positive regulation of neuron projection...
```

Any help would be appreciated?

My entire code:

```
library(topGO)
library(GO.db)
library(biomaRt)
library(Rgraphviz)
library(tidyverse)
library(dbplyr)

prebg_genes <- read.table("mart_mmu_protein_coding_genes.csv", header=F, sep=",")
bg_genes <- as.character(prebg_genes[,1])

pre_candidate_list_1 <- read.table("All_targets_ENSMUSG.txt", header=T, sep=",")
candidate_list_1= as.character(pre_candidate_list_1[,1])

db <- useMart('ENSEMBL_MART_ENSEMBL',dataset='mmusculus_gene_ensembl', host="www.ensembl.org")
go_ids <- getBM(attributes=c('go_id', 'ensembl_gene_id', 'namespace_1003'), filters='ensembl_gene_id', values=bg_genes, mart=db)
go_ids <- read.table("go_ids_mmu_protein_coding_genes.txt", header=T, sep="\t")
listAttributes(db)
View(go_ids)

gene_2_GO <- unstack(go_ids[,c(1,2)])

pre_keep1 <- candidate_list_1 %in% go_ids[,2]
keep1 <- which(pre_keep1==TRUE)
candidate_list_final1 <- candidate_list_1[keep1]

geneList1 <- factor(as.integer(bg_genes %in% candidate_list_final1))
names(geneList1) <- bg_genes

GOdata1 <- new('topGOdata', ontology='BP', allGenes = geneList1, annot = annFUN.gene2GO, gene2GO = gene_2_GO)

classic_fisher_result1 <- runTest(GOdata1, algorithm='classic', statistic='fisher')

weight_fisher_result1 <- runTest(GOdata1, algorithm='weight01', statistic='fisher')

applied to the topGOdata object.
allGO1 <- usedGO(GOdata1)
all_res1 <- GenTable(GOdata1, weightFisher=weight_fisher_result1, orderBy='weightFisher', topNodes=length(allGO1))

p.adj1 <- round(p.adjust(all_res1$weightFisher,method="BH"),digits = 4)

pre_all_res_final1=cbind(all_res1,p.adj1)
all_res_final1 <- pre_all_res_final1[order(pre_all_res_final1$p.adj1),]
write.table(all_res_final1,"summary_topGo_analysis_all_DE_miRs_targets.csv",sep=",",quote=F)
```

Answer: I found the answer in Github. So basically we have to add ```numChar=1000``` to ```all_res1``` i.e:

```
all_res1 <- GenTable(GOdata1, weightFisher=weight_fisher_result1, orderBy='weightFisher', topNodes=length(allGO1), numChar=1000)
```

numChar: truncates GO term descriptions at 1000 chars (basically, disables truncation)
Question: There is a lot of interest in tissue-specific and differentially expressed genes and a lot of people are using RNA-seq to determine if genes are up-regulated in one tissue vs. another. However, I am interested in a particular gene that, based on qPCR, is highly expressed in a particular aphid tissue (bacteriocyte), but it is not differentially expressed (it is also highly expressed overall in the body and in other tissues, like head and gut). So I'm not interested in using RNA-seq to measure differential expression of this gene, because I don't necessarily expect to find it. I am interested in using RNA-seq to confirm high expression in bacteriocyte (relative to other genes in the same gene family). Is there a way to do this? I know that RNA-seq does not measure absolute expression, but would it be valid to rank some measure of expression (like TPM) of, say, my gene of interest and other members of the same gene family within bacteriocyte? I would like to know if this gene and its orthologs are, or are not, highly expressed in bacteriocyte tissues relative to other genes in the same gene family. In theory, I would like to be able to quantitatively compare across orthologs/different insects, but I would be happy with independent yes/no answers for each insect.

Answer: This is really a worst-case scenario for RNAseq. I would suggest that you just do qPCR, which will be both faster and vastly cheaper.

The issue with trying to do this with RNAseq is that there are GC and other biases that will affect the quantification of one isoform versus another. This can still be done, but you really need to be very very comfortable with the data and the underlying biases. It's much easier to do absolute quantification with qPCR, which will directly give you an answer.
Question: Can two MESH ID associated with two different sub type overlap or have same MESH ID .

For example

Acute M2 Myeloid Leukemia	C1879321	M2	D015470

Acute Undifferentiated Leukemia	C0280141	M0	D015456,D015470

This `

> D015470

` is same for both the subtype.

Answer: The answer is yes. In this case, the MeSH term is a more general term for the two more specific concepts. In general, mappings between ontologies are not always perfect 1:1 mappings.
Question: Here is the link to the UniProt ID A0N5G1_HUMAN:
https://www.uniprot.org/uniprot/A0N5G1

I could maximum find this :https://www.ebi.ac.uk/ena/data/view/AAB25742\

Could not find much on this. Any leads will be appreciable.
I am assuming this gene is submitted but under review at transcript level, hence we do not find any associated gene symbol, but I may be wrong. Any insights will be very much appreciated.

P.S.: Posted the link of this Biostar query in Twitter to see if someone responds to the thread.

Answer: Got your query from Twitter :)

This is a TrEMBL [entry][1], meaning an unreviewed sequence with low quality and scored poorly according to UniProt [association score][2] rules. It gives no results in Ensembl (at least for [human][3]), Entrez [Gene][4] or [Open Targets Platform][5], which leads me to believe it's not been used for gene/transcript annotation.

I have also downloaded [all the genes associated with rheumatoid arthritis][6] hoping to get the gene name picked up from text mining or elsewhere but no luck. And [LINK][7] (that allows us to explore connections for genes from millions of PubMed abstract documents and semantic relationships) but still no luck. You may want to delve deeper in the relations for [Rheumatoid][8] only in LINK though.

How about contacting [UniProt][9] helpdesk to clarify on this?

[1]: https://www.uniprot.org/uniprot/A0N5G1#entry_information
[2]: https://www.uniprot.org/help/annotation_score
[3]: https://www.ensembl.org/Human/Search/Results?q=A0N5G1;facet_feature_type=;site=ensembl;facet_species=Human
[4]: https://www.ncbi.nlm.nih.gov/gene/?term=A0N5G1
[5]: https://www.targetvalidation.org/search?src=q:A0N5G1
[6]: https://www.targetvalidation.org/disease/EFO_0000685/associations
[7]: http://blog.opentargets.org/2018/01/18/link/
[8]: https://link.opentargets.io/view?query=rheumatoid
[9]: https://www.uniprot.org/contact
Question: Hi,
I just got a list of genes which were highly expressed in my cells.
I want to determine which genes are transcription factors.
Now I am searching the genes one by one but it needs a long time.
So is there a better and faster way?
Such as a database, where I input all the genes, then the transcription factors will be identified automatically.

Thanks!

Bing

Answer: Yes there are resources to do that. You have to search the google in a much better and fancy way. Data mining is one of the most important feature one should learn. Having said that take a look at this [biostar][1] post to find answer to your query.

You can also se here

1. [TFCat][2]
2. [PAZAR][3]
3. [AnimalTFDB][4]

you can use any of them , to get a list of TFs in Humans if you are working in Humans or other species. But in anycase try to read which was was updated latest to have more robust answer.

[1]: https://www.biostars.org/p/53590/
[2]: http://www.tfcat.ca/
[3]: http://www.pazar.info/
[4]: http://www.bioguo.org/AnimalTFDB/species.php?spe=Homo_sapiens
Question: Hi All,

Does anyone know of a resource in R to obtain the genomic coordinates of cytogenetic bands. I have tried bioMart without success and was hoping to obtain this information through an R interface without resorting to an outside resource.

Answer: You could probably run the following command-line search via a `system()` or MySQL library call within R:

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e "SELECT chrom, chromStart, chromEnd, name FROM cytoBand"

For instance:

> result <- system("mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e \"SELECT chrom, chromStart, chromEnd, name FROM cytoBand\"", intern=TRUE)

Or you might run the RMySQL library:

```
> install.packages("RMySQL", type="source")
> library(RMySQL)
> connection <- dbConnect(MySQL(), user="genome", host="genome-mysql.cse.ucsc.edu", dbname="hg19")
> result <- dbGetQuery(conn=connection, statement="SELECT chrom, chromStart, chromEnd, name FROM cytoBand")
```

The output from `system()` may require additional parsing with `split()` and other functions to build a matrix or data frame.
Question: Does the multiIntersectBed function of the bedtools work with the parameter of fraction (`-f`). I am trying to intersect 4 bed files together to see if they have overlaps together but minimum bases should be around 50% , am not interested in single base match. I am interested in finding if 50% or above of bases in each file match among all bed files giving me the pattern showing the listing of the column where I can see all the labels together. I am trying to do something like this but it does not seem to work. if anyone can give me some suggestions how to get this done? It seems that `-u` and `-f` does not work here. Is there a way I can do that? I tried BEDOPS as well, not that much of a help. I would appreciate some input.

multiIntersectBed -i file1.bed file2.bed file3.bed file4.bed -header -names A B C D -header -u -f .50| head

Regards

Answer: You can't do this directly with bedtools multiinter, but you can answer this by combing a few tools. It is a smidge complicated, but it should work and it demonstrates how combining multiple bedtools can usually get you the answer you need.

```
➜ cat a.bed
chr1 15 20 a1
chr1 18 25 a2
➜ cat b.bed
chr1 16 20 b
➜ cat c.bed
chr1 17 21 c
```

Step 1: Find intervals common to all three files. The fourth column of the output is the number of datasets present at the interval, so it can be used as a filtering criterion.

```
➜ bedtools multiinter -i a.bed b.bed c.bed | awk '$4 == 3' > common.bed
➜ cat common.bed
chr1 17 21 3 1,2,3 1 1 1
```

Step 2. The issue now is that you don't know the degree of overlap these intervals have with the original intervals in each file. So, we will now just use the intersect tool to intersect the common intervals with the original 3 files while requiring 50% overlap.

```
➜ bedtools intersect -a common.bed \
-b a.bed b.bed c.bed \
-f 0.50 -r \
-wa -wb \
> common.with.originals.bed
➜ cat common.with.originals.bed
chr1 17 21 3 1,2,3 1 1 1 1 chr1 15 20 a1
chr1 17 21 3 1,2,3 1 1 1 1 chr1 18 25 a2
chr1 17 21 3 1,2,3 1 1 1 2 chr1 16 25 b
chr1 17 21 3 1,2,3 1 1 1 3 chr1 17 21 c
```

Step 3. Now we have to filter the output from step 2 to ensure that the "common" intervals indeed overlapped with all 3 datasets. Conveniently, when using the `-wb` option with multiple "B" files, the intersect tool reports the dataset from which the hit came. In this case, that will be the 9th column. So, we just need to test that each interval from step 2 had hits from all three input files. We can do this with the groupby tool.

```
➜ bedtools groupby -g 1-8 -c 9 -o count_distinct -i common.with.originals.bed
chr1 17 21 3 1,2,3 1 1 1 3
```

Since the last column now reflect the count of distinct datasets that overlapped the common interval at 50%, you can use awk to filter for only those intervals where the count is 3.

```
➜ bedtools groupby -g 1-8 -c 9 -o count_distinct -i common.with.originals.bed \
| awk '$9==3'
chr1 17 21 3 1,2,3 1 1 1 3
```
Question: I was trying to add a title to an upset plot using the suggestion in this [post][1] but it returned an error.

Error in upset(plot, nintersects = 8, nsets = 6, sets = c("Pos.Corr.K1", :
non-numeric argument to binary operator

[1]: https://www.biostars.org/p/284993/

My code is

library (Upset)
library(grid)
#setting up some random data
Pos.Corr.K1 <- c("GO:0055114", "GO:0016491",
"GO:0016620", "GO:0005739",
"GO:0005759", "GO:0022900",
"GO:0016208")

Neg.Corr.K1 <- c("GO:0005840", "GO:0003987",
"GO:0019427", "GO:0005524",
"GO:0005737", "GO:0043231",
"GO:0000166")

Pos.Corr.DXB11 <- c("GO:0055114", "GO:0016491",
"GO:0005739", "GO:0005737",
"GO:0043231", "GO:0000166",
"GO:0005524")

Neg.Corr.DXB11 <- c("GO:0016620", "GO:0005759",
"GO:0022900", "GO:0005840",
"GO:0003987")

#Draw Upset Plot
data <- list( Pos.Corr.K1, Neg.Corr.K1,
Pos.Corr.DXB11, Neg.Corr.DXB11)
names(data) <- c( "Pos.Corr.K1", "Neg.Corr.K1",
"Pos.Corr.DXB11",
"Neg.Corr.DXB11")
plot <- UpSetR::fromList(data)
setwd("C:/PETER PROJECT/3g. Mass Spec Processing (Comparison of K1 v DXB11)/Results/Upset Plot")
# create PNG for the heat map
png("Upset Plot (2) - K1 v DXB11 (DE Proteins).png",
width = 10*300, # 10 x 300 pixels
height = 10*300,
res = 300,) # 300 pixels per inch
Upset <- upset(plot,
nintersects = 8,
nsets = 6,
sets = c( "Pos.Corr.K1", "Neg.Corr.K1",
"Pos.Corr.DXB11", "Neg.Corr.DXB11"),
#group.by = "sets",
#cutoff = 0,
empty.intersections = "off",
order.by = "freq",
decreasing = T,
point.size = 5,
line.size = 1,
mb.ratio = c(0.65, 0.35),
number.angles = 0,
text.scale = c(2, 2, 2, 2, 2, 2)
#text.scale = c(intersection size title,
# intersection size tick labels, set size title,
# set size tick labels, set names,
# numbers above bars)
) +
grid.text("Distribution of DE Proteins", x = 0.65, y = 0.95,
gp = gpar(fontsize = 20))
print(Upset)

The code runs perfectly without the lines

+
grid.text("Distribution of DE Proteins", x = 0.65, y = 0.95,
gp = gpar(fontsize = 20))

but this is the code suggested in the original post.

Any advice would be appreciated.

Answer: If you run upset() on its own, and run grid.text() afterwards, you get a plot with a title (I do, using your code). Too bad upset doesn't take a title argument - apparently one might be available in UpSetR_1.4.1.

Anyway, you can draw your plot and then you should be able to use dev.copy2pdf() or given your current set up: dev.off() to create the figure.

You also have an extra comma in your png() function call.
Question: I have a bed file as follow:

chr1 100194350 100194710 ARID3A 1000
chr1 151430604 151430964 ARID3A 1000
chr1 20301327 20301687 ARID3A 1000
chr1 229267393 229267573 ARID3A 1000
chr1 8802108 8802375 ARID3A 1000
chr1 109289093 109289349 ATF1 1000
chr1 110527180 110527436 ATF1 1000
chr1 110950342 110950486 ATF1 1000
chr1 115124275 115124409 ATF1 1000
chr1 115259380 115259491 ATF1 1000
...

and I would like to convert it to a GRangeList object in R as follow:

load("data.rda")
head(data)
$ARID3A
GRanges object with 8999 ranges and 0 metadata columns:
seqnames ranges strand

1 chr1 [1307917, 1308277] *
2 chr1 [1407080, 1407440] *
3 chr1 [1858670, 1859030] *
4 chr1 [2175900, 2176260] *
5 chr1 [2290655, 2291015] *
... ... ... ...

8966 chrX [154495495, 154495855] *
8967 chrX [154799333, 154799693] *
8968 chrX [154819952, 154820312] *
8969 chrX [154840885, 154841245] *
8970 chrX [155434904, 155435264] *

-------

seqinfo: 23 sequences from an unspecified genome; no seqlengths

$ATF1

GRanges object with 14883 ranges and 0 metadata columns:
seqnames ranges strand

1 chr1 [ 778593, 778805] *
2 chr1 [1000794, 1001007] *
3 chr1 [1032962, 1033218] *
4 chr1 [1109781, 1110037] *
5 chr1 [1185572, 1185828] *
... ... ... ...

14846 chrX [155026863, 155027119] *
14847 chrX [155057436, 155057692] *
14848 chrX [155881105, 155881361] *
14849 chrX [155881673, 155881929] *
14850 chrX [155893620, 155893876] *
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths

...

I checked GenomicRanges but I could not find a way to make it from bed file?

Thank you so much for helping me.

Answer: We need to split then use lapply to get list output of ranges

# dummy data
df1 <- read.table(text = "chr1 100194350 100194710 ARID3A 1000
chr1 151430604 151430964 ARID3A 1000
chr1 20301327 20301687 ARID3A 1000
chr1 229267393 229267573 ARID3A 1000
chr1 8802108 8802375 ARID3A 1000
chr1 109289093 109289349 ATF1 1000
chr1 110527180 110527436 ATF1 1000
chr1 110950342 110950486 ATF1 1000
chr1 115124275 115124409 ATF1 1000
chr1 115259380 115259491 ATF1 1000 ", header = FALSE)

library(GenomicRanges)

# split and convert per region
res <-
lapply(split(df1, df1$V4), function(i){
GRanges(seqnames = i$V1,
ranges = IRanges(start = i$V2,
end = i$V3,
names = i$V4))
})

# result
res

$ARID3A
GRanges object with 5 ranges and 0 metadata columns:
seqnames ranges strand

ARID3A chr1 [100194350, 100194710] *
ARID3A chr1 [151430604, 151430964] *
ARID3A chr1 [ 20301327, 20301687] *
ARID3A chr1 [229267393, 229267573] *
ARID3A chr1 [ 8802108, 8802375] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths

$ATF1
GRanges object with 5 ranges and 0 metadata columns:
seqnames ranges strand

ATF1 chr1 [109289093, 109289349] *
ATF1 chr1 [110527180, 110527436] *
ATF1 chr1 [110950342, 110950486] *
ATF1 chr1 [115124275, 115124409] *
ATF1 chr1 [115259380, 115259491] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths