cannin/biostars_qa huggingface dataset:
Question: My lab has started using the Haloplex technology to capture the regions of interest.

With this technology most of the reads have the same start/end (I've found that removing the duplicates would reduce the coverage to ~1 base). The reads are grouped in a vertical 'cluster' that can contains more than 1000 reads.

Have you ever used this technology ? I'm currently using samtools mpileup with the option
"-A" ( count anomalous read pairs) and "-d 8000" ( max per-BAM depth).

I'm afraid there is a kind of bias with haloplex: is it safe to use samtools or Gatk to call the variants ?

Pierre

Answer: Agilent will say that you should not de-duplicate Halo data. The vertical clusters are due to the placement of the restriction enzyme sites used to create the capture kits. There's very little off-target capture with Halo and the coverage can vary widely across an capture region ('Manhattan skyline' I think of it as).

You can use, and I have used samtools, GATK or VarScan to call variants from Halo data. However if you do use GATK then if you're post-filtering UnifiedGenotyper calls based on assumptions for exome or WGS data, you will probably find variants which do not pass your standard filtering test that are almost certainly genuine. Strand bias filters etc. seem to be the most prone to being tripped.

I've spoke to Agilent technical support before Christmas about what they propose for best practice genotyping, and the impression I got is that it's something they're working on but don't have a definitive answer for yet. I hope you're trimming the reads in conjunction with their guidelines as well.

EDIT: Since I wrote this post SureCall was released:

http://www.genomics.agilent.com/en/NGS-Data-Analysis-Software/SureCall/?cid=AG-PT-154

This is capable of analysing HaloPlex data and comes directly from Agilent to support the analysis of Haloplex designs from the SureDesign wizard.
Question: Hello all,

I want to calculate the length of sequences from a huge FASTA file. I have read names of which I need to calculate the length from FASTA file. I know how to do it all reads using the following command.

awk '/^>/ {if (seqlen) print seqlen;print;seqlen=0;next} {seqlen+=length($0)}END{print seqlen}' file.fasta

However, I am not sure, how to do it for a specific number of reads. Could someone help me here?

Answer: I would just subset the sequences you need by name using any one of dozens of tools that have been posted on here before, in to a different file, and then just calculate all the lengths in the new file.
Question: I have a hierarchical tree in Newick format, such as:

(A:0.556705,(B:0.251059,C:0.251059):0.305646):0.556705;

I need to collapse clades (sub-trees) whose average distances to the tips (terminal nodes/tips/leaves) are less than a given value, x. Such that input and collapsed output are in Newick format. For example, if B and C in the above collapsed into BC, we get something like:

(A:0.556705,BC:0.556705):0.556705;

Is there a simple way to programmatically achieve this collapsing given any Newick tree (e.g. in Python)?

**Post-Answer Edit**:

Thanks to https://www.biostars.org/u/4370/ and https://www.biostars.org/u/30/ 's answers - for anybody who needs to know whether their solutions actually work or not, they converge to the same outcome given some very complex trees! Bravo

Answer: I had a similar script using ETE. I just wrote a quick example to fit your question:

from ete2 import Tree
def mean(array):
return sum(array)/float(len(array))

def cache_distances(tree):
''' precalculate distances of all nodes to the root'''
node2rootdist = {tree:0}
for node in tree.iter_descendants('preorder'):
node2rootdist[node] = node.dist + node2rootdist[node.up]
return node2rootdist

def collapse(tree, min_dist):
# cache the tip content of each node to reduce the number of times the tree is traversed
node2tips = tree.get_cached_content()
root_distance = cache_distances(tree)

for node in tree.get_descendants('preorder'):
if not node.is_leaf():
avg_distance_to_tips = mean([root_distance[tip]-root_distance[node]
for tip in node2tips[node]])

if avg_distance_to_tips < min_dist:
# do whatever, ete support node annotation, deletion, labeling, etc.

# rename
node.name += ' COLLAPSED avg_d:%g {%s}' %(avg_distance_to_tips,
','.join([tip.name for tip in node2tips[node]]))
# label
node.add_features(collapsed=True)

# set drawing attribute so they look collapsed when displayed with tree.show()
node.img_style['draw_descendants'] = False

# etc...

# Example
t = Tree("((A,(B:0.1,C:0.1)i1:0.1)i2:0.5,((D:0.1,(E:0.1,F:0.1)i3:0.1)i4:0.5,G)i5:0.3);", format=1)
print t.get_ascii(attributes=["dist", 'name'])
# /-1.0, A
# /0.5, i2
# | | /-0.1, B
# | \0.1, i1
# | \-0.1, C
#-1.0, NoName
# | /-0.1, D
# | /0.5, i4
# | | | /-0.1, E
# \0.3, i5 \0.1, i3
# | \-0.1, F
# |
# \-1.0, G

# Now we collapse
collapse(t, 0.5)
print t.get_ascii(attributes=['name'])
# /-A
# /i2 COLLAPSED avg_d:0.466667 {C,A,B}
# | | /-B
# | \i1 COLLAPSED avg_d:0.1 {C,B}
# | \-C
#-NoName
# | /-D
# | /i4 COLLAPSED avg_d:0.166667 {F,D,E}
# | | | /-E
# \i5 \i3 COLLAPSED avg_d:0.1 {F,E}
# | \-F
# |
# \-G

# tree visualization will hide collapsed nodes
t.show()

# collapsed nodes are labeled, so you locate them and prune them
for n in t.search_nodes(collapsed=True):
for ch in n.get_children():
ch.detach()
print t
# /-i2 COLLAPSED avg_d:0.466667 {C,A,B}
#--|
# | /-i4 COLLAPSED avg_d:0.166667 {F,D,E}
# \-|
# \-G

print t.write()
# and the the pruned newick
#(i2 COLLAPSED avg_d_0.466667 {C_A_B}:0.5,(i4 COLLAPSED avg_d_0.166667 {F_D_E}:0.5,G:1)1:0.3);
Question: Dear all,

I have bedg file created from bedtools -

chr1 1 10 num_of_reads

chr1 10 20 num_of_reads

chr1 20 30 num_of_reads

.

.

.

3th column is information about overlapping reads at this region. I would like to add to 4th column information about GC content in reads from my BAM file. Do you have any idea how to do that? Thank you for any help.

So result should be:

ch start stop num_of_reads cg_content

Answer: Here's a python solution. You need to have the [pysam][1] library installed (`pip install pysam`):

```
import pysam
from sys import argv

bam_file = argv[1]
bed_file = argv[2]

bam = pysam.AlignmentFile(bam_file, 'rb')
with open(bed_file) as bf:
for line in bf:
line_parts = line.strip().split()
chr = line_parts[0]
start = int(line_parts[1])
end = int(line_parts[2])
read_data = bam.fetch(chr, start, end)
total_bases = 0
gc_bases = 0
for read in read_data:
seq = read.query_sequence
total_bases += len(seq)
gc_bases += len([x for x in seq if x == 'C' or x == 'G'])
if total_bases == 0:
gc_percent = 'No Reads'
else:
gc_percent = '{0:.2f}%'.format(float(gc_bases)/total_bases * 100)
print '{0}\t{1}'.format(line.strip(), gc_percent)
```

I didn't understand your comment to my question. If you want to print the total number of GC bases in the reads of a coordinate set, just change the last line to:

print '{0}\t{1}'.format(line.strip(), gc_bases)

You run the program like this:

python [script_name] [bam_file] [bed_file]

It prints the modified BED data to STDOUT.

[1]: https://github.com/pysam-developers/pysam
Question: We have a matrix of miRNA expression data (example below w/ random fake data):

```
miRNA sample 1 sample 2 sample 3
hsa-miR-10a-3p 1280 2340 4321
hsa-mir-10-5p 45 48 50
hsa-mir-10b-5p 32 19 25
```

What is the simplest way to generate a heat map + dendrogram based on the jensen-shannon metric using the above data? Do we have to perform clustering and then generate a heat map, or is there an R library that incorporates everything?

EDIT:

Provided an example of what we are looking for below:

Answer: The R heatmap() and heatmap.2() functions can take row and column dendrograms as input. So something like this could work:

clust<-hclust(rowJSdist)
heatmap(data,Rowv=as.dendrogram(clust))

Or you could pass heatmap() the function that computes the distance (and also the one that does the clustering if you need something different from hclust).
Question: I am having a scored bed file with weighted values in the 4th column. I want to use this scored file to as a sorting parameter for another .bam file. However, the input to computeMatrix operation allows only bigWig files to be used in -S option...and creating bigWig files required indexing based on chr# due to which the ranking is lost.

How can I sort the heatmap with a scored BED file ?

computeMatrix reference-point --referencePoint TSS -a 10 -b 10 -bs 2 -p max --sortRegions keep -R UCSC_Main_on_Human_refGene_hg19_TSS.txt -S 8min_sorter_final.bw -out 8min_sorter_final.txt.gz

plotHeatmap -m 8min_sorter_final.txt.gz \
-out 8min_sorter_final.png \
--colorList white,blue \
--missingDataColor 0 \
--heatmapHeight 40 \
--sortRegions no \
--heatmapWidth 10

This gives completely unsorted heatmap.

I am not sure if I am missing something obvious, but I really cant see a way to use a sorted BED file as a direct parameter to sort the matrix generated by computeMatrix. FWIW, I am using Deeptools 2.

Answer: A "scored BED file" is presumably a bedGraph file.

Anyway, what exactly are you trying to plot in the end? If you just want the BAM density over your genes, but the genes sorted according to some sort of score, then you'll need to presort your BED file according to that score. You then don't need the "scored BED file" at all in deepTools. `--sortRegions keep` (this is the default if I remember correctly) in computeMatrix and `--sortRegions no` in plotHeatmap will then result in a plot with that same ordering.
Question: I used following commands to call denovo variants from trio

#Generate a three-sample mpileup
samtools mpileup -B -q 1 \
-f ref.fasta \
dad.bam mom.bam child.bam > trio.mpileup

#Run VarScan trio
varscan trio trio.mpileup trio.mpileup.output \
--min-coverage 10 --min-var-freq 0.20 --p-value 0.05 \
-adj-var-freq 0.05 -adj-p-value 0.15

But, here the output generated from mpileup is really huge. Are there any ways to feed compressed files to the `varscan trio` ?

Answer: You can try: `< zcat (file.gz)` as input

varscan trio < (zcat trio.mileup.gz) trio.mpileup.output \
--min-coverage 10 --min-var-freq 0.20 --p-value 0.05 \
-adj-var-freq 0.05 -adj-p-value 0.15
Question: I know most of you would suggest me to use either `bedtool` or `bcftools`, but I tried a lot with these but could not get full information in my subset vcf file where I intend to extract the SNPs with all column information and updated `AC`, `DP` values etc. in `my_bed.bed` file. When I tried this with `bedtool`, I couldn't process the file due to conflicting AC= column values (it does not change in the subset data), So I am trying to get this done using `vcftools`. However, `vcftools` is also not immune to this problem as I am still missing most of the `FORMAT` column information (`AC=0;AF=0.00054;AN=312;BaseQRankSum=-2.079;DP=40153;Dels=0.00;FS=1.621;HRun=0;HaplotypeScore=0.2550`)

The command I am using to extract using vcftool is :

for file in *.vcf.gz; do
vcftools --gzvcf $file --recode --recode-INFO-all --bed /dir_for_bed_file/my_bed.bed --out "${file}".vcf
done

Can someone please clarify the fuss here or if I am missing anything in my command line?

Answer: Ok figured it out: It was missing `--recode --recode-INFO-all`.
Question: Hi there,

I am working with a VCF that includes Copy Number Variants and hoped to confirm that I am reading the data correctly.

```
22 16050654 DUP_gs_CNV_22_16050654_16063474 A ,,, 100 PASS AC=9,87,599,20;AF=0.00179712,0.0173722,0.119609,0.00399361;AN=5008;CS=DUP_gs;END=16063474;NS=2504;SVTYPE=CNV GT 3|0 0|0 0|0 0|0 0|0
```

From the line above extracted from my VCF I am assuming the following:

- Samples with the 0 genotype have an A at position Chromosome 22 Position 16050654
- For samples exhibiting a CNV the copied block runs from Position 16050654 to 16063474
- Samples with, for example, the 3 genotype have `` copies of this block
- `` means the block is copied 3 times

The final assumption is my main concern, does `` always mean there are x copies of the block?

If anyone can confirm that I am reading this correctly or otherwise set me straight I would be most grateful!

As an aside, as this is my first post, I'd love to thank the Biostars community for all the ongoing help and expertise! You guys are great!

Cheers
MrGraeme

Answer: You should find a definition of CN0,CN1, etc... in the VCF header in the ALT section.

```
Symbolic alternate alleles for imprecise structural variants:
##ALT=
```
Question: I am learning survival analysis in R, especially the Cox proportional hazard model. I read a paper talking about using 80% of the sample as training set and 20% of sample as test set.

As quoted

> On the training set, we first performed a pre-selection step to keep the top significant features correlated with overall survival (univariate Cox model, likelihood ratio test, P< 0.05). ...We used two computational methods to train the models: (i) Cox: the Cox proportional hazards model with LASSO for feature selection ...We then applied the models thereby obtained to the test set for prediction, and calculated the C-index using the R package survcomp.

I do not know how they actually did to apply the models from Cox model to the test set. I mean, for the training set, I can simply perform a coxph function. But the returned results are "coef,exp(coef),se(coef)),z,p" and likelihood ratio test p-value. How can I treat this as a model and use it on the 20% test set data?

Answer: You can try to do something like this:

```
# Derive model in the training data (after feature selection - I believe that in the paper you mentioned they use LASSO: R has a good package for this: glmnet)
cox_model = coxph(Surv(training_data$Survival,training_data$Status) ~ ., data=training_data)

# Create survival estimates on validation data
pred_validation = predict (cox_model, newdata = validation_data)

# Determine concordance
cindex_validation = concordance.index (pred_validation, surv.time = validation_data$Survival,
surv.event=validation_data$Status, method = "noether")
```

See more [here][1]

Hope it helps

[1]: http://stats.stackexchange.com/questions/48298/computing-c-index-for-an-external-validation-of-a-cox-ph-model-with-r
Question: Hi,

This question might be really stupid, but I'm a beginner and information on this would help a lot: I would like to know what exactly people mean when they refer to "steady state mRNA". Could someone help please?

Thanks,
Rama

Answer: from this paper:

> We assume that the change of mRNA level (dR/dt) depends linearly on mRNA synthesis and degradation,

> dR/dt = txj[DNA] - dj[RNAj]

> where [RNAj] is the RNA concentration for gene j, [DNA] is constant ([DNA]=1), txj is the transcription rate, and dj is the degradation rate of gene j. For simplification, we initially assume the degradation rate to be constant, meaning independent of gene j. Therefore in steady state where dR/dt=0, the RNA concentration of gene j is proportional to transcription and degradation rates of gene j.
Question: Hello! everyone:

I'm new to linux, here I got a problem:

I have a file file1 like:

3
6
7
9
12

and file2 which is tab-delimited:

chr1 3052600 3052800 1 E3
chr1 3052800 3053000 2 E3
chr1 3059400 3059600 3 E3
chr1 3059600 3059800 4 E3
chr1 3059800 3060000 5 E3
chr1 3062600 3062800 6 E3
chr1 3101000 3101200 7 E3
chr1 3105000 3105200 8 E3
chr1 3105200 3105400 9 E3
chr1 3116800 3117000 10 E2
chr1 3117000 3117200 11 E2
chr1 3164800 3165000 12 E2

and I want to extract the lines in file2 which its 4-th column equal the number in file1 like below:

chr1 3059400 3059600 3 E3
chr1 3062600 3062800 6 E3
chr1 3101000 3101200 7 E3
chr1 3105200 3105400 9 E3
chr1 3164800 3165000 12 E2

I have spent several hours including wrote a very slow python script, and I searched for the oneline solution, but I got nothing!

awk -v FS="\t" 'NR==FNR{rows[$1]++;next}(substr($NF,1,length($NF)-1) in rows)' fiel1 file2

Thanks a lot for some suggestions!

Answer: This is actually a job for join

join -t $'\t' -1 1 -2 4 -o 2.1,2.2,2.3,2.4,2.5 \
<(sort -g file1) \
<(sort -t $'\t' -k4,4g file2) \
> out
Question: Hey,

I have to do some variant calling in a plant organism (diploid)

I have selected the programs I want to test before I use my own data. (Freebayes, Samtools, SOAPsnp and Cortex_VAR.)
However i am struggling to come up with a reference to compare the tools against.
First i want to test the software before i can use it on my own data.

Is there a benchmark i can use or is it better to simulate my own reference? and use the tools accordantly

Regards,
Covux

Answer: There are a few public reference genomes listed here: http://googlegenomics.readthedocs.io/en/latest/use_cases/discover_public_data/genomic_data_toc.html
Question: Hi,

I am relatively new to bioinformatics/R/scripting. Currently I am trying to improve my coding style so that my scripts are more accessible, both for me and other readers.

Does anyone know of any resources for model scripts I could use as inspiration?

Thanks.

Answer: These are only a few coding style guides I found:

- https://www.bioconductor.org/developers/how-to/coding-style/
- http://r-pkgs.had.co.nz/style.html
- https://web.stanford.edu/class/cs109l/unrestricted/resources/google-style.html
- http://shape-of-code.coding-guidelines.com/2013/03/02/the-most-worthwhile-r-coding-guidelines-i-know/

If you want to look at examples from bioinformatics, you could take a look inside BioConductor packages, like [BioStrings][1]. But these are using advanced features of R's object system like methods, generics, etc. which might be overwhelming for a beginner.

[1]: https://github.com/Bioconductor-mirror/Biostrings/tree/release-3.3
Question: In paired-end RNA-seq analysis, some tools require the insert size (also called inner distance) as the parameter to run (such as [MATS][1]'s `-r` option). And I also know that some tools can be used to estimate the insert size, such as [CollectInsertSizeMetrics][2] in the picard or [inner_distance.py][3] in the RSeQC.

I used both the two tools for my 2*100bp RNA-seq data, however, the results are different, as shown Figure1 in the bellow. It seems the distributions are same, but the values of average insert size are different. It seems that it is because, the CollectInsertSizeMetrics calculates the length of template (the RNA fragment), in contrast, the inner_distance.py calculates the length of template minuse the length of two reads (as shown in Figure2 below, source from RSeQC website).

**Is there anyone know the defination of insert size in RNA-seq? It should inculde the read length (the method of CollectInsertSizeMetrics), or not (the method of inner_distance.py)?**

The mean value of CollectInsertSizeMetrics is 175.825835bp

The mean value of inner_distance.py is -38.1662853494061

(I don't know how to upload a figure into the biostars here.)

Figure1 the mean value and distribution of insert size for my paired-end RNAseq analysis

Figure2 the insert size that inner_distance.py calculates, source for RSeQC website.

[1]: http://rnaseq-mats.sourceforge.net/user_guide.htm
[2]: http://broadinstitute.github.io/picard/command-line-overview.html#CollectInsertSizeMetrics
[3]: http://rseqc.sourceforge.net/#inner-distance-py

Answer: Insert size includes read length

```
----------> <----------
|_______ insert size __________|
```

"Remember that "insert" refers to the DNA fragment between the adaptors, and not the gap between R1 and R2." (http://thegenomefactory.blogspot.de/2013/08/paired-end-read-confusion-library.html)
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: - [Biojava][1]
- [Seqmonk][2]
- [GSEA][3]

[1]: http://biojava.org/wiki/Main_Page
[2]: http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/
[3]: http://software.broadinstitute.org/gsea/index.jsp
Question: Hi,

This appears to be a simple problem that I am unable to solve. I have some data that looks like this:

CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T,G MISSENSE 0.17, 0.09

The above is dummy meaningless data, but it is representative of the problem at hand.

I'd like to `separate_rows` such that the `ALT` and `AF` are separated in a couples manner. Running `separate_rows` on the 2 columns would give me 4 rows, not 2. I'd like my output to be:

CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T MISSENSE 0.17
chr2 1 A G MISSENSE 0.09

Is there any way I can conserve this combination while separating the values out? I am really far out from the VCF to go back and split multi-allelics.

Answer: Using *data.table*:

library(data.table)

x <- fread("CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T,G MISSENSE 0.17,0.09")

x[, lapply(.SD, function(x) unlist(tstrsplit(x, ",", fixed = TRUE))),
by = .(CHROM, POS, REF, TYPE)
][, .(CHROM, POS, REF, ALT, TYPE, AF = as.numeric(AF))]
# CHROM POS REF ALT TYPE AF
# 1: chr1 1 A T MISSENSE 0.23
# 2: chr2 1 A T MISSENSE 0.17
# 3: chr2 1 A G MISSENSE 0.09

Below should work with auto type conversion, but it fails, as the first value "T" gets converted as logical "TRUE", next "T,G" as character, then when binding it, unfortunately, errors out:

x[, lapply(.SD, function(x) unlist(tstrsplit(x, ",", fixed = TRUE, type.convert = TRUE))),
by = .(CHROM, POS, REF, TYPE)]
# Error in `[.data.table`(x, , lapply(.SD, function(x) unlist(tstrsplit(x, :
# Column 1 of result for group 2 is type 'character' but expecting type
# 'logical'. Column types must be consistent for each group.

-------

Related SO post with other alternative solutions:

- [Split comma-separated strings in a column into separate rows](https://stackoverflow.com/q/13773770/680068)
Question: From what I understood, [allele frequency spectrum][1] is everything that concerns the frequency of all of the alleles studied in a population. It can be estimated with Tajima's D, theta, etc.

But how is it different from the *joint* allele frequency spectrum?

[1]: https://en.wikipedia.org/wiki/Allele_frequency_spectrum

Answer: Joint allele frequency spectrum is allele frequency distribution calculated between two\three populations. Imagine a 2D plot, where on x axis is AFS of population 1 and on y axis is AFS in population 2. For example, you can plot derived allele counts in Yoruba versus Central European population. Let's say, you count allele A x times in Yoruba and y times in Central Europe, then you plot it on 2D space. For three populations it will be a 3D plot. Upscaling this to a higher number of populations is computationally very hard.
Joint AFS is informative about migrations and other parameters, because changes in jAFS will reflect changes in demography.

Take a look at this paper by Gutenkunst et al http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1000695