cannin/biostars_qa huggingface dataset:
Question: Hi all,

I am analyzing some Illumina paired-end sequencing experiment. I would like to track the duplicates in my lanes and be able to distinguish between PCR duplicates and optical duplicates.

To this purpose, I use Picard MarkDuplicates.
This function has an OPTICAL_DUPLICATE_PIXEL_DISTANCE parameter ... nice ... but as the function simply set a flag to true in the sorted BAM file, there is no way in the end to distinguish between the two. (Am I right ?)

So, basically I am wondering if this option is really useful ?
It is explained that MarkDuplicates starts to find the 5' coordinates and mapping orientations of each read pair, thus to look at the coordinates of the cluster on the flowcell seems unnecessary (?), as the pair will be tagged as a duplicate anyway.

Do you use in-house script or a particular API for such a goal ?

Cheers
Tony

EDIT : I am aware that Picard creates a metrics file to report some values. But in some lanes generated with a PCR-free protocol, I expected a proportion of my duplicates to be optical duplicates. Nevertheless, in Picard metrics file, I always have %optical_dup=0. So I am wondering if some of you had some issues with this measure as well.

Answer: I omitted to have a close look to the read names in my files. A read name has the following format :

@identifier:lane:tile:x:y

Picard, by default, only match numbers and letters in the 'identifier' part. So if you have underscores (and it's quite usual to have some actually), Picard will not be able to get the coords back and then no optical duplicates will pop up...

Use the READ_NAME_REGEX option of MarkDuplicates to customize the read name matching.
Question: I am new to programming in R and Python, however have some basics. I have a technical question about computation. I would like to know if there are any functions for performing subtraction of all features(rows) to a particular value (row) from the same data list. I would like to obtain the output_value1 as shown below and post this, multiply by (-1) to obtain the output_value2.

Please let me know if you need more details.

I have tried performing the same operation in the MS Excel, this is very tedious and time consuming.

I have many large datasets with several hundred rows and columns which becomes more complex to manually perform the same in MS Excel. Hence, I would prefer to write a code and obtain the desired outputs.

|Feature| |Value| |Output_value1| |Output_value2|

|Gene_1| |14.25633934| |0.80100922| |-0.80100922|

|Gene_2| |16.88394578| |3.42861566| |-3.42861566|

|Gene_3| |16.01| |2.55466988| |-2.55466988|

|Gene_4| |13.82329514| |0.36796502| |-0.36796502|

|Gene_5| |12.96382949| |-0.49150063| |0.49150063|

|Normalizer| |13.45533012| |0| |0|

Answer: Try this example:

# example data
df1 <- read.table(text = "
Feature	Value
Gene_1	14.25633934
Gene_2	16.88394578
Gene_3	16.01
Gene_4	13.82329514
Gene_5	12.96382949
Normalizer	13.45533012
", header = TRUE)

# find the row with normaliser, use the value to substract from each row
df1$Output_value1 <- df1$Value - df1[ df1$Feature == "Normalizer", "Value"]

# then flip the sign (could be skipped, if we flip the substract in above step)
df1$Output_value2 <- df1$Output_value1 * -1
Question: I have a fastq file with over 2 million reads. I am trying to use FastqSampler to select 2 million reads at random. But strangely, I don't get 2 million, I get something less -- 1929951 in the following example. Why?

Might it have something to do with the way FastqSampler chunks the input file? (author is describing this here: A: Selecting random pairs from fastq? )

It works fine if I set n to 1 million reads.

> library(ShortRead)
> fq=readFastq("file.fq")
> fq
class: ShortReadQ
length: 2198402 reads; width: 100 cycles
>
> fqs=FastqSampler("file.fq", n=2e6)
> yield(fqs)
class: ShortReadQ
length: 1929951 reads; width: 100 cycles
>
> sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-redhat-linux-gnu (64-bit)

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

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

other attached packages:
[1] ShortRead_1.14.4 Rsamtools_1.8.6 lattice_0.20-13 Biostrings_2.24.1 GenomicRanges_1.8.13 IRanges_1.14.4

loaded via a namespace (and not attached):
[1] Biobase_2.16.0 grid_2.15.2 hwriter_1.3 tools_2.15.2

Answer: Sorry for the delayed reply. I'm not able to reproduce this with the current version of ShortRead (1.16.3); there were bug fixes related to FastqSampler between the version you are using and the current version.

> set.seed(123)
> fl <- "~/b/working/tmp.fastq"
> readFastq(fl)
class: ShortReadQ
length: 2198402 reads; width: 37 cycles
> fqs=FastqSampler(fl, n=2e6)
> yield(fqs)
class: ShortReadQ
length: 2000000 reads; width: 37 cycles
> sessionInfo()
R version 2.15.2 Patched (2012-12-23 r61401)
Platform: x86_64-unknown-linux-gnu (64-bit)

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

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

other attached packages:
[1] ShortRead_1.16.4 latticeExtra_0.6-24 RColorBrewer_1.0-5
[4] Rsamtools_1.10.2 lattice_0.20-13 Biostrings_2.26.3
[7] GenomicRanges_1.10.7 IRanges_1.16.6 BiocGenerics_0.4.0
[10] BiocInstaller_1.8.3

loaded via a namespace (and not attached):
[1] Biobase_2.18.0 bitops_1.0-5 grid_2.15.2 hwriter_1.3
[5] parallel_2.15.2 stats4_2.15.2 zlibbioc_1.4.0

Try

source("http://bioconductor.org/biocLite.R")
biocLite(character())

to update packages. Please ask follow-up questions on the Bioconductor mailing list (no subscription required).
Question: One thing I really like about CWL is the ability to load CWL files into a CommandLineTool or Workflow object after generating the classes using `schema-salad-tool --codegen=python CommonWorkflowLanguage.yml > cwl_classes.py`. Is there a schema file that describes CWL inputs/job files similar to CommonWorkflowLanguage.yml? I have been working with these files as dicts, but would be super nice to be able to load them straight into a class object.

Answer: Hello Karl,

Yes, the `inputs` section of a CWL document is a schema for the input job object.
Question: Hello !

I want to do a construction and analysis of regulatory genetic networks in specific disease based on involved microRNAs, target genes and transcription factors

Please can you inform me about free tools that can help me to build and analyze a miRNA/Target genes/TFs network.

Thank you in advance.

Answer: We are now building a database and online dataviz tool on cytoscape web -- MirOB (MicroRNA OncoBase).

Database includes validated microRNAs and Transcription Factors Targets associated with human cancers. Currently, MirOB supports only homo sapiens organism.

Network shows all interactions containing transcription factor of miRNA.

You can try our simple tool here.
Question: I am seemingly stuck with something that should be very simple and I hope I haven't overlooked something obvious.

Question: How can I make a valid Blast-database with Taxids from a NCBI query export?

What I have tried so far:

For a meta-genomics project I need a custom made blast database which I wish to generate from the result of the following NCBI Nucleotide query:

Viruses[Organism] AND srcdb_refseq[PROP] NOT cellular organisms [ORGN]

The result is 3986 entries which I exported and saved (via 'Send to') in FASTA and ASN1 format. (Both files are seemingly containing the right amount of entries) As this is a meta-genomics project I would love to have the taxon ids in the blastdb.

I was successful with making a valid blast database from the FASTA file using makeblastdb, but the FASTA header doesn't include taxids, hence I tried to make a blast database from the ASN1 export using the following command (it is not clear from the documentation which formats can be used to create the database):

$ makeblastdb -in AllViralDNARefSeq.asn1 -dbtype nucl -out ViralASN1 -title "All Viral RefSeq DNA from NCBI ASN1"

Building a new DB, current time: 12/20/2011 10:37:28
New DB name: ViralASN1
New DB title: All Viral RefSeq DNA from NCBI ASN1
Sequence type: Nucleotide
Keep Linkouts: T
Keep MBits: T
Maximum file size: 1073741824B
Adding sequences from FASTA; **added 10 sequences** in 0.00906897 seconds.

As you can see, this does not work as it adds only 10 sequences.

Any help to get in the Taxonids is appreciated it doesn't have to be elegant, I just need the database from that query. I am using Blast+ 2.2.25

Answer: I think this one requires some scripting. Here are a few ideas.

First, I'd fetch viral Refseq sequences in Genbank format from the FTP site:

wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.genomic.gbff.gz
gunzip viral.1.genomic.gbff.gz
grep -c "^LOCUS" viral.1.genomic.gbff # 3984

Then I'd parse the Genbank file to extract an ID, description and sequence for each entry. Since taxon IDs are contained in a field of the form `/db_xref="taxon:NNNN` where NNNN = taxon ID, they can be extracted and written into the header of a new fasta file.

Some quick and dirty Bioperl to illustrate:

#!/usr/bin/perl -w

use strict;
use Bio::SeqIO;

my $file = "viral.1.genomic.gbff";
my $seqio = Bio::SeqIO->new(-file => $file, -format => "genbank");
my $fasta = Bio::SeqIO->new(-file => ">$file.fa", -format => "fasta");

while(my $seq = $seqio->next_seq) {
my $taxid = "";
for my $feat($seq->get_SeqFeatures) {
if($feat->has_tag("db_xref")) {
for my $id($feat->get_tag_values("db_xref")) {
if($id =~/taxon:(\d+)/) {
$taxid = $1;
}
}
}
}
my $fa = Bio::Seq->new(-id => $seq->id,
-desc => $taxid.", ".$seq->description,
-seq => $seq->seq);
$fasta->write_seq($fa);
}

This gives a fasta file *viral.1.genomic.gbff.fa*, with headers that look like:

>NC_003038 176652, Invertebrate iridescent virus 6, complete genome.
Question: This question was inspired by Peter Cock 's tweet

https://twitter.com/pjacock/status/118012750105546752

Most softwares (C, java...) use a int32 (unsigned or not) to store the length of the chromosome. It isn't enough when the length of a chromosome is greater than INT_MAX or UINT_MAX

# define INT_MAX 2147483647
# define UINT_MAX 4294967295U

So my question is:

Is there any resource where one can find the length of the chromosomes. Something like BioNumbers.

What's the length of the longest chromosome ?

Answer: I just remembered reading an article on the largest genome size. It claims the flower named Paris Japonica has 150 billion bases over 40 chromosomes,

I would imagine that the largest chromosome is bigger than the UINT_MAX. Alas a quick search did not turn up any sequence information. Most likely the size of the DNA was measured by other means.
Question: let's say I have a fasta of a protein sequence

> albumin

MKWVTFISLL FLFSSAYSRG ... ... ...

I want to split the sequence into all possible consecutive 8 amino acids (only in 1 direction, amino -> carboxyl) (And no looping(I don't know if it is the right expression), say, GMKWVTFIS)

I need

> fasta.albumin1

MKWVTFIS

>fasta.albumin2

KWVTFISL

> fasta.albumin3

WVTFISLLF

...

> fasta.albumin13

FSSAYSRG

And, I want to do this for all known human protein sequences.

How would I do it???
I need the result as a fasta or fasta files. And IDs for resulting 8-mer seuqeunces need to be unique.

Answer: You requirement is to generate simple kmers from sequences.
BioPython solution.
Please change last print statement to: print('>'+str(record.id)+'|kmer_'+str(count)+'\n'+str(my_kmer))

#!/usr/bin/env python3

from Bio import SeqIO

myfile=SeqIO.parse('test.fa','fasta')
for record in myfile:
sequence=record.seq
seq_len=len(sequence)
#define the kmer len
kmer=8
count=0
for seq in list(range(seq_len-(kmer-1))):
count=count+1
my_kmer=(sequence[seq:seq+kmer])
print('>'+str(record.id)+'|kmer_'+str(count)+'\n'+str(my_kmer))
Question: Hello,

I have some microarray data from a few years ago, and the rs identifier in the file was from dbSNP 131. In order to do some analysis on this old data, I will need to know the positions of the microarray SNPs, and thus would like to download dbSNP131.

However, on [UCSC Table Browser][1], I can no long find this database. I tried: assembly: hg19, group: Variations - then there is no dbSNP131 under tracks.

Does anyone know any other places to access dbSNP131 please? Or, is it save to just use dbSNP 138? The microarray chip was HumanOmni2.5-8v1 from Illumina.

Thank you

[1]: http://genome.ucsc.edu/cgi-bin/hgTables

Answer: Okay after some digging, I found that you can downloading all old version of dbSNP databases from [ANNOVAR][1].

Command:

annotate_variation.pl -buildver hg19 -downdb -webfrom annovar snp131 humandb/

The format is exactly the same with the dbSNP table from UCSC: https://genome.ucsc.edu/cgi-bin/hgTables

[1]: http://annovar.openbioinformatics.org/en/latest/
Question: I want to extract **gene name** , **gene start position** and **gene stop position** from the fasta header of the fasta file. I have tried to extract based on the position but those locations are not consistent. Is there any other way to extract them ?

This is what I have tried so far.

#I have a vector of these file names. Here I have just one element

names1 =>"lcl|NC_005336.1_cds_NP_957781.1_1 [locus_tag=ORFVgORF001] [db_xref=GeneID:2947687]
[protein=ORF001 hypothetical protein] [protein_id=NP_957781.1] [location=complement(3162..3611)] [gbkey=CDS]"

#Then I extracted words from the string list

string_list1 <- str_extract_all(names1, boundary("word"))

#result
string_list1[1]

[[1]]
[1] "lcl" "NC_005336.1_cds_NP_957781.1_1"
[3] "locus_tag" "ORFVgORF001"
[5] "db_xref" "GeneID"
[7] "2947687" "protein"
[9] "ORF001" "hypothetical"
[11] "protein" "protein_id"
[13] "NP_957781.1" "location"
[15] "complement" "3162"
[17] "3611" "gbkey"
[19] "CDS"

So, I was trying to extract 4th ,16th and 17th element from this list. It works for this particular example. This does not work for other headers where these positions are different. Usually, gene name is consistently present at the 4th position. But, the start and stop location differ among the fasta headers. So, this strategy is not working and I can't think of any other strategy.

Answer: If gene name is like [locus_tag=gene_name] and coordinates like [location=complement(3162..3611)]

library(tidyverse)

names1 <- "lcl|NC_005336.1_cds_NP_957781.1_1 [locus_tag=ORFVgORF001] [db_xref=GeneID:2947687][protein=ORF001 hypothetical protein] [protein_id=NP_957781.1] [location=complement(3162..3611)] [gbkey=CDS]"

(res <-
str_replace_all(names1,
"^.*?locus_tag=(.*?)\\].*?\\[location.*?(\\d+)\\.\\.(\\d+).*?$",
"\\1___\\2___\\3") %>%
str_split("___")
)
If names will be a "gene_name" column in a data.frame called df, a clean final table can be easily produced:

df %>%
mutate(gene_name = str_replace_all(gene_name,
"^.*?locus_tag=(.*?)\\].*?\\[location.*?(\\d+)\\.\\.(\\d+).*?$",
"\\1___\\2___\\3")) %>%
separate(gene_name,
sep="___",
into = c("gene", "start", "end"))
Question: Hello.

I read the paper "Whole-exome sequencing identifies a recurrent NAB2-STAT6 fusion in solitary fibrous tumors" briefly.

After reading it, I tried to find the gene fusion in my exome sequence data with cancer patients from hospital.

But I think that it is a better strategy to identify the gene fusion case through the use of the public data that is specified in the above paper before trying to test on my own set.

So I tried to find an explanation about how to access public data in the paper to identify the gene fusion.

Finally I found the bold sentence (below).

**Accession code. Sequence data used for this analysis are available at the database of Genotypes and Phenotpyes(dbGap) under accesssion. phs000568.v1.p1.**

Just reading it, I think I can download data if i have permission. So I go to dbGap stite and search phs000568.v1.p1 and many results occur to me. But I don't know how to access above databases and also don't know how to download.

So Is there anyone who have experienced this before?

If so, I need your help!

Thank you.

Answer: http://www.ncbi.nlm.nih.gov/sra?linkname=gap_sra_all&from_uid=874173
Question: Hey all :]

I use samtools's depth, and occasionally samtool's pileup commands to calculate coverage of my reads to the genome before binning for coverage. I'm pretty sure everyone else does too ;)

One very common problem is that people find "samtools depth" and "samtools mpileup" don't match up - and it's commonly attributed to filtering poor quality reads, duplicates, etc (mpileup doing more filtering).

But there's a ton of other questions I just can't get the answers too from the sametools docs, namely:

Does depth/pileup count the region between a paired reads, or just the reads themselves?

Does samtools depth/pileup count singletons in paired end sequencing? Does it extrapolate based on average read length to fillout the whole fragment?

If a read maps to multiple locations, is it counted multiple times?

I'm working with ChIP-Seq data, so i might have to correct for peak-shift. What do you guys think?

Thanks! :)

Answer: > Does depth/pileup count the region between a paired reads, or just the reads themselves?

samtools depth 'just' runs a multi-mpileup and count the number of bases under each base of the reference. smal l'qlen', deletions and base quality below 'baseQ' are ignored. Here is the core code of 'depth':

if (!(b->core.flag&BAM_FUNMAP)) {
if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP;
else if (aux->min_len && bam_cigar2qlen(&b->core, bam1_cigar(b)) < aux->min_len) b->core.flag |= BAM_FUNMAP;

(...)

while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) {
if (pos < beg || pos >= end) continue; // out of range; skip
if (bed && bed_overlap(bed, h->target_name[tid], pos, pos + 1) == 0) continue;
fputs(h->target_name[tid], stdout); printf("\t%d", pos+1);
for (i = 0; i < n; ++i) { // base level filters have to go here
int j, m = 0;
for (j = 0; j < n_plp[i]; ++j) {
const bam_pileup1_t *p = plp[i] + j;
if (p->is_del || p->is_refskip) ++m;
else if (bam1_qual(p->b)[p->qpos] < baseQ) ++m; // low
}
printf("\t%d", n_plp[i] - m); // this the depth to output
}
}
Question: Hi,

I am trying to subsample from a bam file using the samtools view -s command.
This is working when sampling 50% or lower (-s 42.50, 42 being the seed), but anything higher fails (returns an empty file).

He are the exact commands i use

samtools view -s 0.25 -b chr6_all.bam > chr6_25p.sam

Working.

samtools view -s 0.50 -b chr6_all.bam > chr6_50p.sam

Working.

samtools view -s 0.75 -b chr6_all.bam > chr6_75p.sam

Not working.

I also made sure that 49% is working, but 51% is not.
Any ideas, suggestions, or is this an intented mechanic ? There doesn't seem to by any documentation about the subsampling parameter in samtools docfile.

Thanks.

Answer: Subsampling not working for fractions above 50% is a known bug in samtools 0.1.18. (See [[Samtools-help] Randomized Subsampling Bam File / Subsampling above 50%][1].)

The bug was fixed in [March last year][2]; samtools 0.1.19 contains the corrected version.

[1]: http://sourceforge.net/mailarchive/message.php?msg_id=29174919
[2]: https://github.com/samtools/samtools/commit/8cf329904235d82c6c781113c404f298d06bc2b0
Question: Dear all,

I have one problem, it is only one condition for my easy script for counting of GC content per row

awk 'NR>1{n=length($1); gc=gsub("[gcGC]", "", $1); print gc/n}' $i

How to count length of row without N character.

for example:

Input: ACAGCTTGCNNNN => length= 9 Gc content=5/9

format of output is not important, only how to count it.

Thanks a lot

Answer: This command should work, but it's not in awk, it is in `bash`.

```
while read p; do
len=$(echo $p | sed 's/N//g' | tr -d '\n' | wc -c)
cnt=$(echo $p | grep -oh 'C\|G\|g\|c' | tr -d '\n' | wc -c)
gc=$(awk "BEGIN {printf \"%.2f\",${cnt}/${len}}")
echo -e length:$len --- GC:$gc
done
Question: I'm relatively new in this area. Is there instances where published public domain data is reanalysed similar to its original context and reported? Are such papers accepted by well known journals?

As technology and our understanding improves we can generate more information but are there any interest to reanalyse published data? How about ethics? Bioinformatics allows this oppurtunity to analyse or extend existing pool of knowledge but I haven't seen much people out there except those involved in large scale projects/collaborations. Please share your experiences or any papers that are relevant.

Answer: Public data is re-used in new publications in a variety of ways, so I would say this is one of the core features of bioinformatics. Of course, the data should then be used either:

to correct old features, or

to show new characteristics.

I would say the latter has much more chances of being published, except if the corrections are substantial.

Just for fun, and because I think public data is very important for the development of the bioinformatics field, I will quote what one of my previous advisors would say: "Of course most of what we do should be called 'parasitic bioinformatics'. After all, we keep using public data -- generated by other people -- to answer our own questions!"

I could give you an enormous list of papers re-analysing public data :-)
Question: Dear all,

I have downloaded the GRCh38 reference human genome in fasta format from the ftp.ncbi.nlm.nih.gov following the suggestions of this [post][1] and [this][2]. Now I need to generate a fusion genome thus I need the headers right. The headers of the GRCh38 are in this form:

>chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38
>chr2 AC:CM000664.2 gi:568336022 LN:242193529 rl:Chromosome M5:f98db672eb0993dcfdabafe2a882905c AS:GRCh38
>chr3 AC:CM000665.2 gi:568336021 LN:198295559 rl:Chromosome M5:76635a41ea913a405ded820447d067b0 AS:GRCh38
[...]
>chrUn_GL000218v1 AC:GL000218.1 gi:224183305 LN:161147 rl:unplaced M5:1d708b54644c26c7e01c2dad5426d38c AS:GRCh38
>chrEBV AC:AJ507799.2 gi:86261677 LN:171823 rl:decoy M5:6743bd63b3ff2b5b8985d8933c53290a SP:Human_herpesvirus_4 tp:circular

May I ask what is the format of the header? What the individual fields refer to?

Thank you

[1]: http://lh3.github.io/2017/11/13/which-human-reference-genome-to-use
[2]: https://www.biostars.org/p/318811/#334309

Answer: >chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38

AC accession number https://www.ncbi.nlm.nih.gov/nuccore/CM000663.2

gi: ncbi global identfier https://www.ncbi.nlm.nih.gov/nuccore/568336023

LN: length

rl: Assigned-Molecule-Location/Type	(chromosome, mitochondrial...)

M5 : md5 cheksum https://en.wikipedia.org/wiki/Md5sum

AS: assembly