cannin/biostars_qa huggingface dataset:
Question: The Bam format was developed at a time when I/O (disk / network) was not a limiting factor for genomics data analysis. Currently it is. Especially when you have pipelines like GATK that process BAM files in multiple steps and read and write bam files for every step (instead of in memory streaming).

I know raw sequencing data (SOLiD & PacBio) and other large (scientific, financial etc.) datasets are stored in [HDF5](http://www.hdfgroup.org/HDF5/). And there are initiatives like [CRAM](http://www.ebi.ac.uk/ena/about/cram_toolkit) and [reduced bam](https://www.broadinstitute.org/gatk//events/2038/GATKwh0-BP-4-Compression.pdf).

Still I believe most organizations are using regular bam files. Why is this? Is there a roadmap towards development (and adoption!) of improved formats for storing DNA alignment data?

Answer: Firstly, ReducedBAM is still BAM and it has been discontinued. Secondly, not streaming BAM is the problem of GATK, not BAM. A lot of other tools take the advantage of streaming, which SAM/BAM encourages by design. Thirdly, I think CRAM will replace BAM ultimately. CRAM has all the advantages of BAM (except that it is more complex) plus smaller file size and other additional features. With the official SAM/VCF/tabix libraries, htslib in C and htsjdk in Java, seamlessly supporting both, users will be able to use CRAM wherever BAM is applicable.

IMO, the single top reason that other formats (goby, cSRA, adam and biohdf) did not fly is because they are a bit late to the party. They all have their own advantages, but when most data are already in the BAM format and the vast majority of downstream tools accept BAMs only, there is not much room for a new format unless it is substantially better than BAM - not so easy ;-) CRAM wins out for now mainly because it is so close to BAM conceptually that it is easy for them to have the same APIs.

Other minor reasons include... a) [processing BAM was significantly faster][1], though CRAM is catching up now. b) BAM is relatively simple and was fully spec'd from the very beginning. As a result, it has been implemented independently in a variety of programming languages (C, C++ (multiple times), java, lisp, go, D, haskell, javascript, ...). This greatly helps its adoption. c) the following seems a weird reason - SAM/BAM was first implemented in C with admittedly bad but still usable APIs. For a scripting language, binding to a C library is usually easier than binding to a library in C++, Java, etc. Perl, Python and R bindings played an important role in the early adoption of BAM. After all, most analysts work with scripting languages.

[1]: https://sourceforge.net/p/samtools/mailman/message/28989211/
Question: Greetings!

I am doing some analysis in R using sqldf package. The file contains SQL queries of OrthoMCL. It is working fine for the small files but with the larger files (Size of the file is 26 GB) it started to give me error:

> Error in if (n > 0) c(NA_integer_, -n) else integer() :
>
>missing value where TRUE/FALSE needed
>
> Calls: rbind -> rbind -> structure -> .set_row_names
>
> In addition: Warning messages:
>
> 1: In Make.row.names(nmi, ri, ni, nrow):
>
> NAs introduced by coercion to integer range
>
> 2: In nrow + ni : NAs produced by integer overflow

I am unable to understand the rbind error. I looked at other question (https://www.biostars.org/p/98397/) having similar error but I couldn't figure it out. The code for rbind is written as follows:

BestHit2 <- list()
for (i in 1:length(args)){
b <- read.table(args[i],sep="\t",header=T)
BestHit2 <- rbind(BestHit2,b)
}

Any help would be appreciated!

Answer: What is the purpose of using rbind on list element? Probably you need BestHit2 as a data.frame?

> **From rbind/cbind help page**

> Take a sequence of vector, matrix or
> data-frame arguments and combine by columns or rows, respectively.
> These are generic functions with methods for other R classes.
Question: Dear all,

we have been sequencing pediatric tumors by using Illumina technology (normal-tumor pairs, paired end reads of 150nt each), and I have been using several SV callers (including DELLY, MANTA, LUMPY) in order to call the Structural Variants.

Considering the initial files of SV calls, what is the set of criteria you would recommend for filtering ? particularly, the number of PR (paired-end reads) or SR (split - reads) ? thank you,

-- bogdan

Answer: This depends on the relative cost of false positives vs false negatives. What FDR and sensitivity are you aiming for?

I've done some fairly extensive benchmarking of SV callers for which you can find germ-line and simulation results at http://shiny.wehi.edu.au/cameron.d/sv_benchmark to given you some idea of the trade-off between FDR and sensitivity for various callers. The [DREAM Somatic Mutation Calling Challenge results][1] also give you a basic idea of these trade-offs.

In my somatic usage of GRIDSS, I've found that requiring assembly support has been a very effective filter - something that you could also do for manta. As for direct read support, the relative strength of RP and SR support depends on the library fragment size distribution and will be different for each library. For 2x150bp sequencing, you'll have weak RP support even for real variants in libraries with mean fragment size <300bp. Ideally, you would filter on variant quality directly, instead of RP and SR independently but not many SV callers actually report a variant quality score, and none that I know of report an even remotely calibrated quality score.

Disclaimer: I am the author of [GRIDSS][2].

[1]: http://dreamchallenges.org/project/icgc-tcga-dream-somatic-mutation-calling-challenge/
[2]: http://biorxiv.org/content/early/2017/02/21/110387
Question: Hello all,

I am in a fix. For my work, I need an exhaustive and up-to-date database of

all protein sequences, especially covering ALL eukaryotes sequenced till date.

Thus, I downloaded the latest version of nr database from ncbi's ftp site.

However, I find that it does not contain all the putative protein

sequences as listed in individual genome databases.

For discussion purposes, let us consider Cyanophora paradoxa (taxonomic

id: 2762). According to the website

http://cyanophora.rutgers.edu/cyanophora/blast.php, it has 32,167 protein

coding sequences. However, there are only 731 gi ids corresponding to this

species in the latest nr database (25May2014 version). The Cyanophora

paradoxa's complete genome was published in 2012 (Price DC et al, 2012).

Thus, to me the only option seems to be to download protein sequences from individual genome projects and combining identical entries in one entry. Finally I shall append them to nr database to get the exhaustive set of proteins. However before I begin on this mammoth task, I thought to enquire if there is simpler solution to this problem.

Any help or suggestions are greatly welcomed.

Thanks a lot for your time,

Ramya

Answer: Well the current NCBI nr contains 40,337,612 sequences originating from: GenBank CDS translations (excluding those from environmental samples and WGS projects), UniProtKB/SwissProt, PDB and PRF.

UniProt's UniParc database has more coverage, currently containing 63,875,797 unique protein sequences, including the CDS translations from the INSDC databases (DDBJ, EMBL-Bank and GenBank) which are excluded from NCBI nr, and additional sequences from various other sources.

You could also try having a look at SIMAP, this contains additional protein sequences from meta-genomics experiments.

However looking at the paper (Cyanophora paradoxa genome elucidates origin of photosynthesis in algae and plants) and searching the nucleotide databases it appears that the only submission has been to the Sequence Read Archive:

- http://www.ncbi.nlm.nih.gov/Traces/sra/?study=SRP009206
- http://www.ebi.ac.uk/ena/data/view/SRP009206&display=html

The corresponding assembly and associated feature annotations do not appear to have been submitted as yet, presumably due to this currently being a draft genome. Thus the protein sequences which are derived through translation of INSDC CDS features do not appear in NCBI nr or UniProtKB.

In general you should be able to start from a non-identical sequence archive such as UniParc and add unique protein sequences from other sources, assuming you can obtain them with appropriate annotations for your purposes.
Question: I am doing a DGE analysis of a total RNAseq dataset of 2 timepoints (5 reps each). I am particularly interested in looking for changes in expression of 1000 genes.

Currently, I have done the analysis by analysing all genes and then picking the 1000 genes I am interested in.
However, my PI has suggested that I could try doing the the analysis by doing DGE on just the 1000 genes. Theoretically, this should improve the statistical significance since there would be minimal adjustments for multiple hypothesis testing.

Is this an advisable way of doing the analysis? Since differential expression levels are fit to a negative binomial distribution (in the case of DESeq2), wouldn't this just mean most of the 1000 genes I input would end up not being differentially expressed?

Edit: We arrived at the list of 1000 genes as we were interested particularly in genes coding small proteins. Hence, we searched Uniprot for human proteins with a maximum length of 100 amino acids.

Answer: In my opinion this is not an advisable way of doing the analysis. The main problem is how one arrives at the list of interest. In your case, it seems that these 1000 genes were selected *a posteriori* by their statistical significance and not "biological" reasons. So for me it is hardly justifiable.
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: About R graphics, an introduction by Paul Murrell is very popular. And, you know what? He shares all the figures code and 3 chapters via his website. It's a very good one.
Question: Hi all,

I will be using Drop-seq to prepare cDNA libraries for thousands of single cells, followed by single-end sequencing on a HiSeq2500 (read length = 100 bases). This will involve the addition of a unique 12 nucleotide cell barcode to the 5' end of all reads originating from the same cell (so bases 1-12 of each read). Unique molecular identifiers will also be used (they will be bases 13-20 of the reads). I won't be able to use the pipeline designed by the creators of drop-seq because they require paired-end sequencing. My problem is with demultiplexing: I am using randomly generated cell barcodes which will be supplied in excess, so I have no way of knowing beforehand which barcodes were actually used. Because of this, I am unable to supply the barcode sequence information that most scripts out there require as input. As someone with minimal computational/bioinformatics skills, I am not comfortable writing custom scripts to fit my needs.

Does anyone know of any scripts/packages that I would be able to use to demultiplex my data based on the in-line cell barcodes? I am planning on using the python package UMI tools (https://github.com/CGATOxford/UMI-tools) to extract the UMIs and deduplicate reads, but I have not been able to find any information on how to separate the reads from different cells into distinct fastq files.

Thanks in advance for any recommendations!

Answer: I believe `fastq-multx` can do help. It works with in-read barcodes. The `-n` flag should "print likely barcode list". See: https://github.com/ExpressionAnalysis/ea-utils/blob/wiki/FastqMultx.md

If that doesn't help, you can determine which barcodes you have and extract the barcode sequences yourself (read FASTQ file, get the lines with the nucleotide sequence, take the first 12 characters, sort, keep only unique sequences):

zcat file.fastq.gz | sed -n '2~4p' | cut -c 1-12 | sort | uniq

But that will be noisy. To get the most occurring barcodes (more likely to be real) which you can use for demultiplexing afterwards:

zcat file.fastq.gz | sed -n '2~4p' | cut -c 1-12 | sort | uniq -c | sort -nr | head -100
Question: What is the best / standard place to get a full list of SNPs and their coordinates in hg38?

I downloaded the SNPsnap database, but just realized that those coordinates are in hg19.

I'm trying to figure out how many SNP sites exist in my targeted genome sequencing data.

Many thanks.

Answer: One can download it in many formats by first going here and then choosing the dbSNP build version and the human genome reference build:

- ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/

For example, human_9606_b151_GRCh38p7 is dbSNP release version 151 with co-ordinates for GRCh38.p7.

The VCF format is a common download. In the VCF directory, the *00-All.vcf.gz* file is the one that contains all records. Take a look at the READMEs in order to see what's in all of the other files.

Kevin
Question: Is there a way to tabix a text file and use it for quick retrieval of records?

I have these fields in a text file

#CHR START STOP GENE

I would like to retrieve the records based on chr:start-pos.

My sample file [sample.txt] has this data

```
chr1 11873 14361 A
chr1 14469 29370 B
chr1 34610 36081 C
chr1 69090 70008 D
chr1 134772 140566 E
chr1 323891 328581 F
chr1 367658 368597 G
```

I bgzipped it to get sample.txt.gz file

I tabixed with the following command

tabix -b 1 -e 3 -S 0 -0 sample.txt.gz

and this is the error I encountered

```
[get_intv] the following line cannot be parsed and skipped: chr1 11873 14361 A
[ti_index_core] the indexes overlap or are out of bounds
```

Any ideas on how to tabix a text file?

Answer: 'start' is in column 2 , you used `-b 1`

You'd better use the predefined setting `tabix -p bed`
Question: Hello,

I am trying to write a dictionary object to a FASTA file however I have problems with writing it.
I could not achieve doing it without using the library or with the library (Biopython).

I tried converting my dictionary to list using "dict.items()" then writing it with SeqIOand error is:

"AttributeError: 'tuple' object has no attribute 'id'"

I would appreciate any kind of help.
Thanks in advance!

WM

Answer: The error is clear. Your dictionary is missing the `id` attribute, which is the required parameter to use `SeqIO.write`. Typically, you'd provide a `SeqRecord` object to it, which includes a `Seq` object with parameters like `id` and `description`. You can definitely turn your dictionary into a list of `SeqRecord` objects and everything else should work.

There are quite a few other ways to convert dictionary to the many formats `SeqIO` supports. The easiest (and the least programming experience required) is to simply write your dictionary into a tab-delimited file and use `SeqIO.convert`.

See below for an example.

from Bio import SeqIO

a = {'myseq1':'acgt', 'myseq2': 'gctc'}

# try writing your own code to turn this dictionary into a tab-delimited file (seq.tab), i.e
# myseq1 acgt
# myseq2 gctc

SeqIO.convert('seq.tab', 'tab', 'seq.fa', 'fasta')
Question: Hi,

I have a couple of bed files with regions (peaks of ChIPseq), and I would like to create a unique file with merge regions, where overlapping regions are fused IF the overlap cover more than xx % (in my case 60%) of one of the 2 regions (not necessary both).

I know how to do that with 2 files, however I have more (at least 4). I could do recursive and pipes to do it pairwise multiple times (with bedops or bedtools), but I was wondering if there would not be a cleaner way. I can use bedops to merge everything that overlap by 1bp for my bed files at once, but I do not find a way to control on the overlap fraction...

Any idea?

Answer: Let's look at the following [BEDOPS][1]-based solution, which is a modification of the approach authored [here][2] by Shane Neph:

```
$ bedops --everything file1.bed file2.bed ... fileN.bed \
| bedmap --echo-map --fraction-both 0.6 - \
| awk '(split($0, a, ";") > 1)' - \
| sed 's/\;/\n/g' - \
| sort-bed - \
| uniq - \
> answer.bed
```

It might look a bit overwhelming, but it is really quite straightforward if we walk through this line by line.

1. Take the multiset union of all elements from all *N* input files with `bedops --everything`. Note that this is *not* the same as a `--merge` operation, which instead makes contiguous regions from all overlapping input regions. Pipe this unioned set to `bedmap`.
2. For each element from all files, use `bedmap --echo-map` to print any other elements from the multiset union that overlap that first element by 60% or more, using `--fraction-both 0.6`. The result will be a semi-colon delimited string that contains one or more overlapping elements. Note that this result will always include the original element, because it overlaps itself by 100%. Pipe this result to the `awk` statement.
3. The `awk` command splits the string on the semi-colon delimiter and prints out the string if there are two or more overlapping elements. We do this filter step, because we are not interested in any element overlapping itself and only itself, but only multiple (2+) elements that meet the 60% or greater overlap criterion. We pipe this string to `sed`.
4. The `sed` statement replaces all semi-colons with newline characters. This creates BED-formatted output that we pipe to `sort-bed`.
5. The `sort-bed` tool takes the overlapping elements and puts them in lexicographically-sorted order. We pipe this to `uniq` to remove duplicate elements.

The file `answer.bed` is a sorted BED file that contains all unique elements from `file1.bed` through `fileN.bed`, where there is mutual 60% or greater overlap.

Note that this is a highly efficient solution, compared with O(2^N) pairwise comparisons required by other approaches. Further, we make use of Unix standard input and output streaming (note the hyphen at the end of each command, denoting standard input) to avoid the cost of making unnecessary intermediate files; writing to disk is expensive. The only thing written to disk is the final result.

[1]: http://bedops.readthedocs.org/en/latest/
[2]: http://bedops.readthedocs.org/en/latest/content/usage-examples/multiple-inputs.html
Question: Hi

I am trying to create heatmaps. I have used heatmap.2, pheatmap and clustvis. I have dataset of 2000 rows containing allele frquency of different SNPs in different populations. I cannot filter rows as I am visualizing relationship between different populations on the basis of allele frequencies. The problems with R packages is that they does not create good resolution heatmaps. And clustvis plots heatmaps on the basis of principal components not on the basis of actual allele freqiencies.

Can anyone please tell me how to great good resolution heatmaps of such large datasets?

Answer: The resolution or overall sense of 'quality' of a figure is determined by a few different factors, including (but not limited to):

- your monitor's type and capabilities (LED, etc)
- the screen resolution you have set via your operating system
- the compression format (or not) in which you're encoding your figure
(PNG, TIFF, PDF, etc)
- how much information you are attempting to 'pack' into the figure

Problems with heatmaps in R are usually related to too much information being 'packed' into too small a space.

The first thing that you should do is save your figure in a vector-based format, such as PDF, and then you can manage the 'packing' of information into this by altering the width and height, for example:

pdf("Test.pdf", width=7, height=7)
...
dev.off()

Increase size:

pdf("Test.pdf", width=12, height=18)
...
dev.off()

`width` and `height` here relate to inches.

With other formats, like TIFF, PNG, *et cetera*, width and height relate to pixels, but the resolution will never have the clarity as in vector-based formats.

Kevin
Question: Hello, I am using the locuszoom standalone to generate images for publication.

I need to be able to print them at 300 DPI.

I have looked through [the documentation][1], but cannot seem to find an option to do so, even though this must be possible because such images are frequently published.

Can anyone offer help or advice on how to do this?

[1]: http://genome.sph.umich.edu/wiki/LocusZoom_Standalone

Answer: If your images are in pdf you can use the *convert* tool in [imageMagick][1] like:

convert -density 300 in.pdf out.jpg

For the record, this is the command I happened to use to convert pdf figures to tiff suitable for submission to PlOS One:

convert -alpha off -trim -resize 1180x1180 -density 300 -depth 8 -compress lzw fig.pdf fig.tiff

Options explained:

-alpha off : Turn off transparency of background (Plos doesn't like it)
-trim : Remove unused/white space around figures.
-resize : Resize to i x j pixels. If the source image is not squared, make the
largest dim (h or w) these many pixels and
the other dim smaller. I.e. proportions are conserved.
-density : Density as pixels/inch (DPI).
-depth : Default 16, Plos wants 8
-compress : Compression algorithms. lzw, required by Plos, makes the image much smaller!

[1]: http://www.imagemagick.org/script/convert.php
Question: I am looking for a library to parse a GTF file in C++. I'm using htslib for SAM/BAM, and it works well. However, I am not aware of any library that can handle GTF format for me. Is there a good library that I can use to parse GTF file, or should implement a parser myself?

Answer: I have C code for parsing GTF in *convert2bed* on my Github page: https://github.com/alexpreynolds/convert2bed

The relevant function parses an input buffer into a `c2b_gtf_t` struct. You could leave out the last line of the conversion function, if you do not need BED output.