cannin/biostars_qa huggingface dataset:
Question: I was looking into the samtools and bcftools variant callers. When I ran this

samtools mpileup -C 50 -E -t SP DP -u -I -f -b >

It told me that flags t, u and I were deprecated. So I thought that bcftools wast the new way to use the "mpileup" but I'm not sure if they are the same.

bcftools mpileup --threads 16 -C 50 -E -f -b >

In addition, I realized that the multithreaded flag in the bcftools if *only* to compress the file. The variant caller will still be *single core*.

But my question here is what are the differences between `samtools mpileup` and `bcftools mpileup`?

Answer: `Bcftools mpileup` should be used instead of `samtools mpileup` for variant calling. That is, the VCF / BCF output mode of mpileup is better in bcftools.

`Samtools mpileup` however has two different formats with the default always being a simple columnar format showing chr, pos, reference, depth, base-calls and qualities. It's a simple parseable format used by some tools that has no attempt at variant calling.

Samtools has also supported for a long time VCF / BCF output (eg `samtools mpileup -vu in.bam`) which traditionally was the input to bcftools call. However this has been deprecated since samtools 1.9 (it's not even in the help message, although the option still works) and since then it now produces a warning message (although still works) if you attempt to use it to generate VCF:

```
[warning] samtools mpileup option `v` is functional, but deprecated. Please switch to using bcftools mpileup in future.
```

With the next release it'll probably be removed completely, leaving only the original basic text form.
Question: I have some ``ATAC-seq PE`` samples and `Dnase-seq SE` samples, I like to do peak calling for them using MACS2. I run MACS2 command with `--nomodel -- shift -100 -- extsize 200 ` For DNase samples but for ATAc seq samples I am not sure which one of `--nomodel -- shift -100 -- extsize 200 ` or only use ` -f BAMPE` options of MACS2 is better?

Also, I found in ENDCODE https://docs.google.com/document/d/1f0Cm4vRyDQDu0bMehHD7P7KOMxTOP-HiNoIvL1VcBt8/edit part 2a, that use bed file but there 2 files BEDPEand TAG file

-tagAlign file ${FINAL_TA_FILE}

-BEDPE file (with read pairs on each line) ${FINAL_BEDPE_FILE}

-Subsampled tagAlign file for CC analysis ${SUBSAMPLED_TA_FILE}

that I am not sure which one should be use as Input of MACS2. should I use TAG file as Input?

Answer: The [official Harvard guidelines][1] as of this year suggest the following:

> Our previous recommendation was to run MACS2 with -f BAMPE, which is
> similar to the default analysis mode of Genrich (inferring full
> fragments, rather than cut site intervals). Others have attempted to
> interpret cut site intervals with MACS2 by using the --shift and
> --extsize arguments, but these arguments are ignored in BAMPE mode. They do work in the default (BAM) mode, but then, with paired-end
> reads, most of the alignments are automatically discarded (half of the
> properly paired alignments and all of the unpaired alignments;
> secondary alignments are never considered). Is it worse to interpret
> full fragments that may be less informative biologically, or to
> disregard more than half of the sequence data? A complicated question.
> The correct answer is: use Genrich.

The main issue with `--nomodel -- shift -100 -- extsize 200 -f BAM` is that you lose a lot of reads (approx 50%). As the MACS2 manual states specifically that for the `-f BAM` parameter:

> If the BAM file is generated for paired-end data, MACS will only keep
> the left mate(5' end) tag."

So you lose all your 3' information.

If you are interested you could read the [old Harvard guidelines which do cover the optimal parameters for using MACS2][2].

[1]: https://informatics.fas.harvard.edu/atac-seq-guidelines.html
[2]: https://informatics.fas.harvard.edu/atac-seq-guidelines-old-version.html
Question: Can someone verify if this R code for converting raw counts to TPM is correct?

#' @title Compute TPM for a read count matrix
#' @param dfr A numeric data.frame of read counts with samples (columns) and genes (rows).
#' @param len A vector of gene cds length equal to number of rows of dfr.
#'
r_tpm <- function(dfr,len)
{
dfr1 <- sweep(dfr,MARGIN=1,(len/10^4),`/`)
scf <- colSums(dfr1)/(10^6)
return(sweep(dfr1,2,scf,`/`))
}

Answer: Use [this code snippet][1] from Michael Love.

[1]: https://support.bioconductor.org/p/91218/#91256
Question: A colleague of mine and I have just been discussing the meaning of fold changes and though this question has been asked before non of the answers are actually as straight forward as needed so lets try here to solve this once and for all.

I list here what I understand so far and will update it from your answers.

What is the correct way to understand a fold change value of a gene or protein?

- A foldchange describes the difference of two values (eg. difference of expression in gene/protein A between healthy and diseased case)
- Biostatistical porgrams/packages calculate it via: "Log(FC)" = mean(log2(Group1)) - mean(log2(Group2))
- log2 fold changes are used/plotted in graphs as those are nicer to show because they center around 0, giving reductions a negative value and increments a positive value
- log2 fold change values (eg 1 or 2 or 3) can be converted to fold changes by taking 2^1 or 2^2 or 2^3 = 1 or 4 or 8
- To convert the fold change into change in % or anything that is actually tangible/understandable in "real life terms" ... need answers here! (= actual question I want to ask)

Lets get this solved once and for all, im looking forwards to your posts!

Here are great posts explaining more about fold changes:

[graphical representations][1]

[conversion of log2 fold change to fold change][2]

[1]: https://www.biostars.org/p/312980/
[2]: https://www.biostars.org/p/100460/

Answer: There is what appears to be an error in your definition of fold change (could be an error, could be causal wording):

* A fold change describes the **ratio** of two values (not the difference). i.e. (expression condition 1)/(expression condition 2)
* The log2 fold changes are the log-of-the-fold-changes i.e. log2(condition1/condition2)
* Because log(A/B) = log (A) - log(B), many statistical programs will calculate the Log2FC = log2(condition1) - log2(condition2), but this is mathematically identical to Log2FC = log2(condition1/condition2)
* Biostats programs will often estimate log2(condition1) using mean(log2(condition1)). This is equivalent to taking the geometric mean of the original data. Thus, Log2FC = mean(log2(condition1)) - mean(log2(condition2)) is the same as Log2FC = log2(geo_mean(condition1)/geo_mean(condition2)).
* Log2 fold changes are used/plotted in graphs as those are nicer to show because they center around 0, giving reductions a negative value and increments a positive value
* log2 fold change values (eg 1 or 2 or 3) can be converted to fold changes by taking 2^1 or 2^2 or 2^3 = 1 or 4 or 8
* You can interpret fold changes as follows:
* **if there is a two fold increase (fold change=2, Log2FC=1) between A
and B, then A is twice as big as B (or A is 200% of B).**
* **If there is a two fold decrease
(fold change = 0.5, Log2FC= -1) between A and B, then A is half as
big as B (or B is twice as big as A, or A is 50% of B).**

@kristoffer.vittingseerup has the formula for converting foldchange or LogFC into % change.
Question: After running vcf-stats on a vcf file I get terms like:

'snp_count',
'private',
'hom_AA_count',
'hom_RR_count',
'ref_count',
'count',
'unphased',
'het_RA_count',
'ref',
'indel_count'

`snp_count` should be the number of snps, and `indel_count` is the number of indels. Private from googling appears to be number of SNPs (or maybe variants) exclusive to a sample in a multi-sample SNP. However, what do the rest of them mean? What is the difference between `ref` and `ref_count`?

Answer: From a few months ago: https://www.biostars.org/p/83558/
Question: I installed ABySS following the instructions from the official documentation.

>sudo apt-get install abyss

and then proceeded with a simple test

>wget http://www.bcgsc.ca/platform/bioinfo/software/abyss/releases/1.3.4/test-data.tar.gz
>tar xzvf test-data.tar.gz
>abyss-pe k=25 name=test in='test-data/reads1.fastq test-data/reads2.fastq'

Everything looks fine except when I run

>abyss-fac test-unitigs.fa
abyss-fac: command not found

It seems that the only avaiable command is abyss-pe.

Any idea of what could be the problem?

Answer: If you look at [the Debian package page][1], the only thing installed to `/usr/bin` is `abyss-pe`, which is why you can find it from the command line. Try:

/usr/lib/abyss/abyss-fac test-unitigs.fa

[1]: https://packages.debian.org/wheezy/amd64/abyss/filelist
Question: Good day. I am new to genome annotation. I am running maker on a new genome and am planning on using Braker to train augustus and gene mark. I have already used my tophat alignment of my mRNA-seq data to the genome to run braker. My question is how then do I use the braker output in maker? Is there a maker script that converts the gtf files produced by braker into a format accepted by maker in the CTL files? Thank you. Will appreciate any help

Answer: BRAKER1 uses GeneMark-ET to train and evaluate AUGUSTUS, and the final gene predictions are from AUGUSTUS.

You can do two things: First, you can use the model that BRAKER creates in the AUGUSTUS config dir in MAKER's AUGUSTUS run. All you have to do is to specify the species name you gave to BRAKER in MAKER's maker_opts.ctl. For example, you ran BRAKER like this:

braker.pl --genome=genome.fasta --bam=reads.bam --species=my_species

This will create a folder named `my_species` in the config dir of AUGUSTUS. If you use the same AUGUSTUS in MAKER all you have to do is enter your species in MAKER's maker_opts.ctl, and AUGUSTUS will find it again:

augustus_species=my_species #Augustus gene prediction species model

Alternatively, you can treat the BRAKER annotation as "legacy annotation" in MAKER. In that case, you can use the final .gff from BRAKER1 as `pred_gff` in MAKER's maker_opts.ctl, and MAKER will edit the exons and introns based on the external evidence you give it. There is an example for this in [this protocol][1].

I have no idea which of the two approaches is "better". If you need more help with running MAKER/BRAKER, the Supplementary Materials from the BRAKER publication has many useful commands: http://bioinformatics.oxfordjournals.org/content/suppl/2015/11/09/btv661.DC1/supplementary.pdf

Later edit: One caveat you will run into is that MAKER performs repeat masking while BRAKER does not, so you're (especially with plant species) bound to see more false-positive transposon-related genes with BRAKER if you don't run repeat masking first.

[1]: http://www.yandell-lab.org/publications/pdf/maker_current_protocols.pdf
Question: I run roary ( roary -e --mafft -p 32 *.gff) to produce core genome alignment on hundreds of Salmonella sequences and have results. I have the following files in the results:

1. gene_presence_absence.csv
2. gene_presence_absence.Rtab
3. pan_genome_reference.fa
4. accessory_binary_genes.fa.newick
5. accessory_graph.dot
6. core_accessory_graph.dot
7. core_gene_alignment.aln
8. clustered_proteins

I would like to visualise the panSNP tree with the roary matrix of these sequences using the command: roary_plots.py name_of_your_newick_tree_file.tre gene_presence_absence.csv. But I realise I have no .tre file in my roary outcome. How do I generate this .tre file from my sequences or from the roary outcome I have?

Answer: you can use a program like FastTree where you input the core_gene.aln file and it outputs the Newick tree file

fasttree core_gene.aln > tree.file

http://www.microbesonline.org/fasttree/
Question: Can anyone tell me how to generate `vcf.gz` file and its index file `vcf.gz.tbi` in 1000 Genomes Project?

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20110521/

Answer: bgzip -c file.vcf > file.vcf.gz
tabix -p vcf file.vcf.gz

[tabix documentation][1]

[1]: http://samtools.sourceforge.net/tabix.shtml
Question: Hi

I am trying to do permutation test for comparing number of LoF mutations in cancer patients.

Haven't done this before based on gene size and GC content. Do you have any statistical tutorial for this? or code?

Thanks

Answer: In R you can use the package [regioneR][1] available in Bioconductor. It takes care of the randomization process and the statistics and you can customize it to evaluate different features.

If I understand correctly your question, you want to test if the proportion of LoF variants you have in you sample is higher that what you would expect by chance, is'nt it? If so, you could try with something like this: create a number of random variant sets based on yours, evaluate how many LoF variant you detect, see if the original set has more LoF variants than expected.

I would do something like this:

1 - Implement a function counting the number of LoF variants given a GRanges with the reference and alternative alleles. You can use the predictCoding function in [VariantAnnotation][2], for example, with this signature:

```r
numberOfLoF <- function(A, ...) {
varAllele <- DNAStringSet(as.character(A$alt))
mcols(A) <- NULL
cod <- predictCoding(query=A, subject=txdb, seqSource=Hsapiens, varAllele=varAllele)

#Note1: the same variant can be counted more than once if it affects more than one transcript
#Note2: What you consider a LoF can be different than that
numLoF <- length(which(cod$CONSEQUENCE %in% c("nonsynonymous", "frameshift", "nonsense")))

return(numLoF)
}
```

2 - Implement a randomization that places the variants randomly over the genome and assigns the same nucleotide changes you see in your variants:

```r
randomizeVariants <- function(A, ...) {
randA <- randomizeRegions(A, ...)
mcols(randA) <- mcols(A)
return(randA)
}
```

3 - Create a mask to tell the randomization function where the random variants can be placed. If this is from an exome sequencing experiment, I would mask out any part of the genome not targeted by the exome.

```r
exome.targets <- toGRanges("path/to/exome.bed")
whole.genome <- filterChromosomes(getGenome(genome="hg19"), chr.type="canonical")
mask <- subtractRegions(whole.genome, exome.targets)
```

4 - And finally, perform a permutation test with all of this

```r
library(regioneR)
library(VariantAnnotation)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

#Define here the functions and mask above

#Permutation test

pt <- permTest(A=vars, randomize.function=randomizeVariants, evaluate.function=numberOfLoF,
mask=mask, genome="hg19", ntimes=50)
print(pt)
plot(pt)
```

This may take some time to run, depending on the number of variants you have and the ntimes, but it will perform the whole permutation test and give you a p-value, a z-score and a plot of the results of the permutation tests.

You can find more information about regioneR in the package [vignette][3].

[1]: http://bioconductor.org/packages/release/bioc/html/regioneR.html
[2]: http://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html
[3]: http://www.bioconductor.org/packages/release/bioc/vignettes/regioneR/inst/doc/regioneR.pdf
Question: Hello!

Recently, I'm learning about RNA-seq quatification and I found there are many normalization methods, such as, CMP,RPKM/FPKM,TPM,TMM etc.
After reading some papers about them, I noticed that all of these units mentioned above can't compare directly between samples because they just get relative expression instead of absolute expression.

I have just heard about spike-in method which may not be a proper method to get absolute expression, so, I want to ask you whether there exsits or you have imaged a proper expression unit to get absolute expression or compare between samples.

Thanks for your suggestion!

Answer: > compare between samples.

Do you want to do differential expression analysis? If so, the most appropriate tools are DESeq2 and edgeR. These will also deal with normalisation and provide the most sophisticated method so you don't have to worry about absolute vs relative quantification.
Question: Hi Biostars Users,

I would like to know if you have encountered this issue while working with [CONTRA tool][1] for detecting CNVs in targeted NGS.

[1]: http://sourceforge.net/apps/mediawiki/contra-cnv/index.php?title=CONTRA:_Copy_Number_Analysis_for_Targeted_Resequencing#Examples

In fact, my analysis is "stopping" at the binning process with no error message, I have checked all the required arguments and they seem to be ok :

```
contra apple$ python contra.py --target ngs1_probes.bed --test cic04668.bam --control baseline_ngs1.bed --bed --fasta human_ref.fasta --outFolder ~/contrases

target : ngs1_probes.bed
test : cic04668.bam
control : baseline_ngs1.bed
fasta : human_ref.fasta
outfolder : /Users/apple/contrases
numBin : [20]
minreaddepth : 10
minNBases : 10
sam : False
pval : 0.05
sampleName : No-SampleName
nomultimapped : False
plot : False
bedInput : True
minExon : 2000
largeDeletion : False
Creating Output Folder : Done.
Converting TEST Sample...
DEBUG 123 genomeCoverageBed -ibam cic04668.bam -bga -g /Users/apple/contrases/buf/sample.Genome
Converting CONTROL Sample...
Getting targeted regions DOC...
chr1
chr10
chr11
chr12
chr14
chr15
chr16
chr17
chr19
chr2
chr20
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrX
Targeted regions pre-processing: Done
Getting targeted regions DOC...
chr1
chr10
chr11
chr12
chr14
chr15
chr16
chr17
chr19
chr2
chr20
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrX
Targeted regions pre-processing: Done
Test file read depth = 97323185
Control file read depth = 98289477
Pre-processing Completed.
Getting the Log Ratio
Binning....
```

Do you have any idea why?

Thank you for your help

Kiz

Answer: Hello Kizauna,

I am new to CONTRA and I am facing the same problem. I have no idea what is wrong as there is no error message. I am sure the thread is not stopped but it is running in my server for a couple of days. Were you able to resolve this.

Thanks in advance
Question: Hello, biostars!
I have some problems with miltyPhylo object manipulation and can't understand how to solve it.

I have many trees from MrBayes analysis and want to use it for comparative phylogenetics methods. Trees are ultrametric and UNrooted (also can't understand why, in MrBayes my outgroup sequence was the first in the matrix).
Than I did following steps:
1) outgroup deletion (species from other group, which I don't want to include in analysis).

>lapply (ntrees [1], is.ultrametric)
$` TREE * gen.6720000`
[1] TRUE
>lapply (ntrees [1], is.rooted)
$` TREE * gen.6720000`
[1] TRUE

Than I looked at one of the trees, it was unrooted.

2) To root trees for the following analysis

>rtrees<-lapply (ntrees, root, "new outgroup")
> lapply (rtrees [1], is.rooted)
$` TREE * gen.6720000`
[1] FALSE

> lapply (rtrees [1], is.ultrametric)
$` TREE * gen.6720000`
[1] FALSE

I'm new in phylogenetic comparative methods and in R project. So I just can't find a mistake.

Answer: Hi Alice,

A couple of things here.

1) the apply family functions are for applying a given function to every member of a list/matrix/vector etc. But you've sliced out the first part of the multiPhyloobject (basically a list). So you are only applying the function to the first part. On top of that, there are two ways to slice lists in R, single or double square brackets. Using single brackets will return a list (with the slice indices you give it), double brackets the particular member of this list you are denoting. (see ?Extract to get your head around the distinction). The has some unexpected results:

trees <- replicate(10, rtree(10), simplify=FALSE)
plot(trees[[1]]) #a tree, works fine
plot(trees[1]) # a list, R doesn't know how to plot it
plot(trees[1][[1]])# the first element of the list is a tree, plots fine

In your example you are applylng the is.rooted and is.ultrametric functions to lists not trees, and that may have similar unexpected results.

To check all trees in a tree object are rooted just drop the slice:

lapply(ntrees, is.rooted) #screeds of TRUEs
all(lapply(ntrees, is.rooted)) #one TRUE

2) You probably need to specify resolve.root=TRUE when rooting the tree with a single tip. This will add 0-length branch to node below the specified tip (which becomes the root). Check out ?root for the details. You can add additional arguments to the function call to apply family functions:

rtrees <- lapply (ntrees, root, "new outgroup", resolve.root=TRUE)
Question: Hello,

I'm new here and and I need help please!

I need to download a list of all human genes with their respective gene symbol | chromosome | strand | transcription start site | Txen | and **Ensembl gene name**,

Actually I'm using UCSC table and I get something like this:

```
#hg38.knownGene.name hg38.knownGene.chrom hg38.knownGene.strand hg38.knownGene.txStart hg38.knownGene.txEnd hg38.kgXref.geneSymbol
uc001aaa.3 chr1 + 11873 14409 DDX11L1
```

And I want to know if it's possible to include ensembl gene symbol with UCSC table or with another method

Thank you in advance

Cherif

Answer: ### [Introduction to BioMart][1]

[Mining data][2]

[Data Mining in Ensembl][3]

[1]: http://www.ensembl.org/Multi/Help/Movie?db=core;id=189
[2]: http://www.ensembl.org/info/website/tutorials/biomart_worked_ex.pdf
[3]: http://www.ensembl.org/info/website/tutorials/mining_with_biomart.pdf
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: I have found these two helpful:

S Poetry by Patrick J Burns

An introduction to S and S plus by Phil Spector