cannin/biostars_qa huggingface dataset:
Question: I'm trying to understand the formats of different SNP callers. I have made test bam files with a subset of two stickleback WGS samples. I made a vcf using VarScan's commad `mpileup2snp`.

Here's an example of a single SNP in the vcf;

chrI 1789 . T G . PASS ADP=105;WT=1;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/0:66:108:108:87:18:17.14%:1.7211E-6:63:74:64:23:15:3 0/1:68:102:102:73:21:22.34%:1.3568E-7:60:56:48:25:19:2

I can understand most of the scores in the `#FORMAT` field, however I'm confused how DP is calculated. The reference and alternative allele depths are given (RD & AD), but the sum of the two doesn't equal DP. Using the example above;

Sample 1:

- DP = 108
- AD = 18
- RD = 87
- AD + RD = 105

Sample 2 :

- DP = 102
- AD = 21
- RD = 73
- AD + RD = 94

Does anyone have any ideas why the reported DP doesn't equal AD + RD?

Answer: Different variant have different rules regarding thresholds for counting reads. Additionally, these criteria may vary for different fields. VarScan offers an [explanation on their website][1]:

> VarScan requires that bases meet the minimum Phred quality score
> (default 15 for most commands) to count them for things like read
> counts (reads1, reads2) and to compute variant allele frequency.
> However, when VarScan reports the depth (such as in the DP field of
> VCF output), it reports SAMtools raw depth. To get VarScan read counts
> to more closely match another tool, set use parameter --min-avg-qual
> 0. And use caution! Low-quality bases, with the occasional exception of BAQ penalties, should not be trusted.
>
> Also, VarScan reports variants on a biallelic basis. That is, for a
> given SNP call, the "reads1" column is the number of
> reference-supporting reads (RD), and the "reads2" column is the number
> of variant-supporting reads (AD). There may be additional reads at
> that position showing other bases (SNP or indel variants). If these
> other variants meet the calling criteria, they will be reported in
> their own line. If not, it may look like you have "missing" reads.

[1]: http://varscan.sourceforge.net/support-faq.html
Question: Hi, this question is somehow related to this previous question. I'm playing with the C API for BAM and I wrote the following code: https://gist.github.com/736059

I'm now trying to find the indexes of the genomic bases covered by the CIGAR string (my final goal is to create a WIG file containing the coverage of the genome)

What would be the correct way to change my code to get the genomic indexes of the bases covered by each CIGAR element ?

(...)
for( k=0;k< b->core.n_cigar;++k)
{
int cop =cigar[k] & BAM_CIGAR_MASK; // operation
int cl = cigar[k] >> BAM_CIGAR_SHIFT; // length
switch(cop)
{
case BAM_CMATCH: printf("M");break;
case BAM_CINS: printf("I");break;
case BAM_CDEL: printf("D");break;
case BAM_CREF_SKIP: printf("N"); break;
case BAM_CSOFT_CLIP: printf("S");break;
case BAM_CHARD_CLIP: printf("R");break;
case BAM_CPAD: printf("P");break;
default:printf("?");break;
}
printf("%d",cl);
}
(...)

Thanks

Answer: For example

```
int depth[REF_LEN], x, j, k;
uint32_t *cigar = bam1_cigar(b);
for (k = 0, x = b->core.pos; k < b->core.n_cigar; ++k) {
int op = cigar[k]&16;
int l = cigar[k]>>4;
if (op == BAM_CMATCH) {
for (j = x; j < x + l; ++j) ++depth[j];
x += l;
} else if (op == BAM_CREF_SKIP || op == BAM_CDEL) x += l;
}
```

Something like this...
Question: Hi,

Could provide the faster way to filter this data :

chr1	43	1000	gene_name=boby	gene_type=trucA
chr2	44	1000	gene_name=natt	gene_type=trucB
chr3	45	1000	gene_name=alurika	gene_type=trucC

To :

chr1 43 1000 boby trucA
chr1 44 1000 natt trucB
chr1 45 1000 alurika trucC

CORRECTION :
Original text data looks like this :

chr1 43 1000 TEST gene_name=boby;gene_type=trucA;foo=34
chr2 44 1000 TRUC gene_name=natt;gene_type=trucB;foo=34
chr3 45 1000 PASS gene_name=alurika;gene_type=trucC;foo=34

Answer: sed -e 's/gene_name=//g' -e 's/gene_type=//g' file > file2
Question: I would like to "mimic" the online BLASTn wgs search (NOT nt, this works already) on my local machine. The reason is that I want to extend from the BLAST nt search to a broader range of bacterial strains.

It seems like there is no BLAST wgs database to download. How do I determine which genomes I have to download to catch the genomes in wgs. I only need it for bacteria for a start.

An alternative would be to not use the entire wgs but some method to at least "broaden" the search compared to BLAST nt.

Answer: Did you take a look at [this README][1] file NCBI provides for instructions on how to deal with WGS data with `blast+`? Scripts you need are in [this directory][2]. Bacteria are `taxid 2`.

[1]: ftp://ftp.ncbi.nih.gov/blast/WGS_TOOLS/README_BLASTWGS.txt
[2]: ftp://ftp.ncbi.nih.gov/blast/WGS_TOOLS/
Question: Hi Biostars Leaders,

Freebayes(version:v1.0.1-1-g683b3cc-dirty) defines AF as Description="Estimated allele frequency in the range (0,1]", but the values are always either 0.5 or 1.0, and they are not actual observed frequencies.

I have observed the same with GATK's HaplotypeCaller, and I have custom calculated the actual frequencies from Ref & Alt Alleles i.e. from the AD field.

Freebayes does not spit out the AD field, but it has these following fields which I think I can use :
RO = "Reference allele observation count, with partial observations recorded fractionally"
AO = "Alternate allele observations, with partial observations recorded fractionally"

I am wondering if there is any advice on how to calculate actual allele frequencies for Freebayes ?

thanks,
gsr

Answer: I got answer to my own question.

I just installed a newer version of freebayes (v1.1.0-3-g961e5f3-dirty) which does spit out the AD field "Number of observation for each allele"

I then used the AD field to calculate the actual allele frequencies in human sample(NA12878)
Question: I have few VCF Files which contains Chromosome number, position, Variation allele (Ref and Alt).

How can I generate the Gene name, its Refseq and listing all of its transcripts?
1.	Any tool which gives these kind of results when upload VCF data?
2.	Is there any module in Biopython or Bioconductor or script for this?
3.	Is there any lists that links gene name to Refseq and transcript information?

Thank you for the valuable suggestions

Answer: You need the [variant effect predictor][1], does exactly what you need.

[1]: http://www.ensembl.org/info/docs/tools/vep/index.html
Question: For drug [Fluvastatin][1], [Fluvastatin [sodium salt]][2] (CAS-93957-55-2) is usually used for testing the potential of Fluvastatin. They have different structures, though it's slight difference.

In [Drugbank][3], the Fluvastatin structure is at the beginning, while the salt form in the latter section of the webpage.

It seems the active form should be Fluvastatin ion in cell.

So which structure should be used when considering the structure ot this drug in chemoinformatics? Thank you.

Update:

Another example:

For [**tamoxifen**][4], one of its Synonyms is tamoxifen citrate shown in guidetopharmacology. In pubchem, they are different (CID [2733526][5] and [2733525][6]) with the latter includes a citrate. At least in the Connectivity Map, tamoxifen citrate [54965-24-1] (as a catalog name) was used instead of tamoxifen (as a cmap_name) [[ref1][7]]. I think they are just different, resulting in different effect on cells. In this situation, should I use the structure of tamoxifen to discover something in the Connectivity Map? It's a little odd. I'm wondering why they use tamoxifen citrate instead of tamoxifen, but name it as tamoxifen. Is it common for drug testing in chemistry?

[1]: https://pubchem.ncbi.nlm.nih.gov/compound/1548972#section=Top
[2]: https://pubchem.ncbi.nlm.nih.gov/compound/23663976#section=Top
[3]: http://www.drugbank.ca/drugs/DB01095
[4]: http://www.guidetopharmacology.org/GRAC/LigandDisplayForward?ligandId=1016
[5]: https://pubchem.ncbi.nlm.nih.gov/compound/2733526#section=Top
[6]: https://pubchem.ncbi.nlm.nih.gov/compound/2733525#section=Top
[7]: http://www.connectivitymap.org/cmap/cmap_instances_02.xls

Answer: You have chosen a tricky one because it has salts, racemates, enantiomers and virtual deuteration

this adds up to 113 different representations http://www.ncbi.nlm.nih.gov/pccompound?cmd=Link&LinkName=pccompound_pccompound_parent_connectivity_pulldown&from_uid=1548972

Choices depends on exactly what cheminformatic operations you intend to perform, but stripping salts back to parents is usually advisable

To keep life simple just go with this GtoPdb entry

http://www.guidetopharmacology.org/GRAC/LigandDisplayForward?ligandId=2951
Question: where to get the file of hg19 exon, intron, UTR region ?

I read lots of post in biostar. However it's been a long time post.

Answer: using **bioalcidaejdk** : http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

```
$ wget -q -O - "ftp://ftp.ensemblgenomes.org/pub/release-45/metazoa/gtf/apis_mellifera/Apis_mellifera.Amel_4.5.45.gtf.gz" |\
gunzip -c |\
java -jar dist/bioalcidaejdk.jar -F GTF -f biostar.code

(...)
6	4682136	4682216	+	GB52198-RA.Intron16
6	4682387	4682473	+	GB52198-RA.Intron17
6	4682696	4682760	+	GB52198-RA.Intron18
6	4682905	4682967	+	GB52198-RA.Intron19
6	4676837	4677042	+	5' UTR of GB52198-RA
6	4683076	4683853	+	3' UTR of GB52198-RA
6	4691339	4691384	+	GB52199-RA.Exon1
6	4692448	4692491	+	GB52199-RA.Exon2
6	4693914	4694249	+	GB52199-RA.Exon3
6	4691384	4692448	+	GB52199-RA.Intron1
6	4692491	4693914	+	GB52199-RA.Intron2
6	4691339	4691339	+	5' UTR of GB52199-RA
(...)
```

with biostar.code:

```
stream().
flatMap(GENE->GENE.getTranscripts().stream()).
flatMap(TRANSCRIPT->{
final List L = new ArrayList<>();
TRANSCRIPT.getExons().stream().forEach(E->L.add(E.toInterval()));
TRANSCRIPT.getIntrons().stream().forEach(I->L.add(I.toInterval()));
TRANSCRIPT.getUTRs().stream().forEach(U->L.add(U.toInterval()));
return L.stream();
}).forEach(R->println(R.getContig()+"\t"+(R.getStart()-1)+"\t"+R.getEnd()+"\t"+R.getStrand()+"\t"+R.getName()));
```
Question: I am trying to pipe the output from BWA to sambamba to sort and index the sam files. I have 20 files with reads from sequencing (pair end) and want to have the resulting bam file (not the intermediate sam or bam files). This is the code I have at the minute:

for filename in ./seqtk_1/subsample_1/*_1.fq.gz;
do file=`echo $filename|sed 's/_1.fq.gz//'`;
filenopath=`basename $file`;
outputpath=./BWA/seqtk_1/subsample_1;
bwa mem -v 3 ./combine_reference.fa.gz ${file}_1.fq.gz ${file}_2.fq.gz > ${outputpath}/align_${filenopath}_BWA.sam |
sambamba view -S -f bam - > ${outputpath}/align_${filenopath}_BWA.bam |
sambamba sort -o - > ${outputpath}/sorted_${filenopath}_BWA.bam |
sambamba index - > {outputpath}/indexed_${filenopath}_BWA.bam;
done

This is the output:

-bash: {outputpath}/indexed_sub_NC_001539_BWA.bam: No such file or directory
sambamba-view: Unrecognized option -
sambamba-sort: Cannot open or create file '' : No such file or directory
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 100000 sequences (10000000 bp)...
[M::process] read 100000 sequences (10000000 bp)...

That continues through the rest of the files. I get a sam file and a sorted_${filenopath}_BWA.bam file but the bam file isnt populated.

My thinking is that the code isn't read/completed linearly and it is trying to create files that can't be created because BWA hasn't started running yet.

Is there a way to fix this? Or do I just need to run BWA and sambamba separately? I don't want to keep these sam files because the size is too large.

Thanks in advance

Answer: I don't know sambamba, so I don't know if you can use it to read from stdin, but I'll assume it does.

Your concept of pipes is wrong. Normally you would do the following

command1 input > output1
command2 output1 > output2

With pipes you would do

command1 input | command2 > output2

So you don't send the output of bwa to a sam file, but directly to sambamba.
Question: In short, my question is:
1. what is the running time of building FM-Index, is it linear in the sense of reference genome?

I was told that FM-Index can be built in linear time, i.e. O(n), where n is the length of reference genome.
But I actually can not find the paper describing it.

The paper cited mostly "Indexing Compressed Text" by Ferragina et al. 2004 focus on the analysis of memory footprint and query time analysis, but not the running time of building FM Index.

Answer: In 2003, two (or maybe three) conference papers first showed that suffix array can be constructed in linear time. As generating BWT from suffix array takes linear time, they also proved that FM-index can be constructed in linear time. One of the most influential linear-time algorithms, SA-IS, was invented by Nong et al (2008). They provided a ~100-line implementation in C. It is much simpler and faster than the previous works. Yuta Mori optimized Nong et al's implementation into the sais library. A few years ago (when I was still following the literature), that library was fastest linear-time implementation in practical benchmarks. However, at least at that time, libdivsufsort, an O(nlogn) algorithm developed also by Yuta Mori, was widely believed to be the fastest open-source library to construct suffix arrays, faster than sais. Linear-time algorithms are not necessarily faster. If a linear algorithm is associated with a large constant and lots of cache misses, an O(nlogn) algorithm can be faster in practice, just as Rob said.
Question: Previously, I split out a vcf file by chromosome, and for my project, I have combined the X and XY vcf files into a single one. After changing the "XY" chromosome designation to "X" via:

awk '{gsub(/"XY"/, "X"); print;}' Genome_newX.vcf > Genome_newX2.vcf

I'm running into the issue of sorting this new "Genome_newX2.vcf" by position. The idea is that I'll subsequently run the vcf through GenotypeHarmonizer.

Are there any suggestions on how to do this easily? I'm brand new to this style of work, and I'd love some direction on where to read up on it as well. Thank you!

Answer: Would there be an equivalent for a BCF? `bcftools view | [...]` code? Or why not using `bcftools sort -Oz output.bcf -o output_sort.vcf.gz`?
Question: After running picard ValidateSamFile I get errors for all reads like the one below - "NM" tags are missing.

WARNING: Read name SRR6251016.24364087_TGTTATGAGA, A record is missing a read group
WARNING: Record 1, Read name SRR6251016.24364087_TGTTATGAGA, NM tag (nucleotide differences) is missing

I am using bam files produced by a STAR mapping pipeline which have "nM" tags as shown below. These are identical in function to NM tags but are alternatively named.

SRR6251016.24364087_TGTTATGAGA 99 chr1 3043025 255 70M = 3043191 236 AGAAAATTGGACATAGTACTACCGGAGGATCCAGCAATACCTCTCCTGGGCATATATCCAGAAGATGCCC EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE<

Answer: If STAR's `nM` was actually identical in function to the standard `NM`, one might ask why on earth they were making life hard for everyone by using a different tag name. But in fact it is not:

> nM : is the number of mismatches per (paired) alignment, not to be confused with NM, which is the number of mismatches in each mate.

Look at STAR's `--outSAMattributes` option, which can be used to also output `NM`. One might ask the STAR developers why `NM` and other tags desired by Picard's “typical usage” validation are not in STAR's standard set of attributes…
Question: Hi,

Assume a table as below:

X =

col1 col2 col3
row1 "A" "0" "1"
row2 "B" "2" "NA"
row3 "C" "1" "2"

I select combinations of two rows, using the code below:

pair <- apply(X, 2, combn, m=2)

This returns a matrix of the form:

pair =

[,1] [,2] [,3]
[1,] "A" "0" "1"
[2,] "B" "2" NA
[3,] "A" "0" "1"
[4,] "C" "1" "2"
[5,] "B" "2" NA
[6,] "C" "1" "2"

I wish to iterate over pair, taking two rows at a time, i.e. first isolate `[1,]` and `[2,]`, then `[3,]` and `[4,]` and finally, `[5,]` and `[6,]`. These rows will then be passed as arguments to regression models, i.e. `lm(Y ~ row[i]*row[j])`.

I am dealing with a large dataset. Can anybody advise how to iterate over a matrix two rows at a time, assign those rows to variables and pass as arguments to a function?

Thanks,
S ;-)

**Edit:** In response to the comments, I should specify that my problem concerns SNP and expression data where I aim to do a pairwise multiple regression analysis (first order regression) in order to assess any possible SNP-SNP interactions that may effect the expression phenotype.

Answer: Using `mapply` you can map a function over more than one sequence (list or vector) in one pass e.g. the odd and even rows of a matrix

fn <- function(i, j) paste("i =", i, "j =", j)

odds <- seq(1, 9, by = 2)
evens <- odds + 1

mapply(fn, odds, evens)
Question: I have two genotype tables, one for parents generated from snp arrays following hg18 coordinates, and the other for their offspring generated from ngs following hg19 coordinates. I finished liftover from hg18 to hg19 for parental data (refer to this post) and I got a new list of coordinates, but I found I couldn't merge the parental genotype with the new hg19 coordinates because I don't know how to pair the new coordinates with the old ones. I wonder if there is a way to get their matching relationship. Any idea? Thanks.

Answer: Thanks, Ashutosh. Problem solved. I changed the 4th column to unique IDs, eg 1.rs#, 2.rs#, 3.rs#. Then my worries are gone.
Question: Dear All,

I wonder if R stops support for "makeTxDbFromUCSC" function? I tried many times but keep failing with the error message below:

**> txDB <- makeTxDbFromUCSC(genome = "hg19", tablename = "ensGene")**

**Error in makeTxDbFromUCSC(genome = "hg19", tablename = "ensGene") :
could not find function "makeTxDbFromUCSC"**

Any input is appreciated.

Thanks,

Xiao

Answer: Did you load the library?

library(GenomicFeatures)

See possible solutions at SO:

- [Error: could not find function … in R](https://stackoverflow.com/q/7027288/680068)
Question: Hi all,

I'm currently working with some GWAS data after imputation using IMPUTE, however, we are having problem to conduct SNPTEST, and has been suggested to convert the GEN and Sample files after IMPUTE to PLINK BED/PED files, then run association analysis using PLINK, to double check.

Using GTOOL, has been able to convert GEN and sample files to PED and MAP files, but when run logistic regression, ERROR: Locus rs has >2 alleles:
individual 1 has genotype [ A G ]
but we've already seen [ N ] and [ A ], occurred. Have tried to remove the relevant individual and snp, but still getting the same message.

The code used for the above message in PLINK was --file test_chr01 --logistic --....--out..

Trying to figure out why before to disturb my supervisor again. Any advices are appreciated. Many thanks.

Guan

Answer: I had the same problem before and one of my group members helped me out.

The error is related to that missing data: GTOOL makes the PED file with 'N' as the missing genotype, which plink doesn't recognize as missing because plink expects missing genotypes to show '-9'. Plink can understand the 'N' for missing using the following command-line parameter:

--missing-genotype N

Add this when converting GEN and sample files to PED/MAP.
Question: Hi,

I need phyloP scpres for the whole human genome. I was able to download 10M positions from http://genome.ucsc.edu/cgi-bin/hgTables using these instructions:

group: Comparative Genomics

track: Conservation

table: your phyloP table of choice

region: genome

output format: data points

However I need the whole genome information, not just 10M.

The website directs me to the downloads page, but I was lost there unfortunately.

Can anyone please link me to the correct place? (Human Grch37 phyloP values in data point format)

Thanks!!

Answer: https://www.biostars.org/p/86847/#86892

https://www.biostars.org/p/242484/#242486
Question: I have gone through filtering variants through `VQSR` and `hard-filter` from [here][1]

My understanding about `VQSR` is that the we don't want to combine `SNP` and `INDEL` where as they are combined in `Hard-filter`.

1. In VQSR, we run `VariantRecalibrator` with mode `SNP` and `INDEL` and we get `.recal` files for both snp and indels.

2. Next, we apply `ApplyVQSR` with mode `INDEL` with `indels.recal` file to generate `indel.recalibrated.vcf`.

3. In the next step, we apply `ApplyVQSR` for `SNP` with vcf input as `indel.recalibrated.vcf` and .recal file generated from `VariantRecalibrator` with mode `SNP`.

4. This step generate file with `snp.recalibrated.vcf.gz` which will
contain both `SNP` and `INDEL` as filtered and will be final filtered data from `VQSR`

Is my understanding correct about variant filtering here ?

If this is right how do we deal with `MIXED` type ?

[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360035531112

Answer: It is better to split the multi-allelic site to bi-allelic sites: use `bcftools`

bcftools norm \
-m-any \
--check-ref -w -f /path/to/reference/hg38.fasta \
input.vcf -o output.vcf
Question: I often have ChIP-seq experiments where I want to get a feel for read coverage over the predicted binding regions. Looking through the UCSC genome browser region by region is slow and labourious. What i am looking for is a way of plotting read coverage (from bedGraph, Wig etc) for many different binding regions and present many of these plots on one page.

I envisage the input would be a bedGraph/Wig file and a list of binding region coordinates. I am aware of a previous Biostars thread that kind of covers this, but it uses UCSC/IGV. I am looking for something much more simplistic - just a line graph per region. Even more, it would be great to be able to plot ChIP and input read coverage on the same graph.

I wonder whether some Python guru etc. has already tackled this?

Many thanks!

Answer: I will like to update my previous answer:

We have made available a suite of tools called deepTools to make this sort of visualizations very easy. There is a tool called profiler that will plot exactly what you want. You provide a list of regions (BED or GFF format) and a bigWig file and the output is a profile. Here you can find the documentation on how to run the tool. You can easily convert any wig/bedgraph file to a bigWig using the tools from UCSC that can dowloaded here: http://hgdownload.cse.ucsc.edu/admin/exe/

The tool uses bigWig files to parallelize computations.
Question: Hello,

Is it possible to parse the read entries in an (unaligned) BAM file one by one, without requiring huge amounts of memory?

I've tried `samtools view | ` as well as using the `pysam.AlignmentFile()` parser from inside the script, but both solutions on our cluster end up taking over 60GB of RAM for a 6GB BAM. I do believe we have nodes that have can handle a lot more ram than that, but I'm still annoyed by requirements that wouldn't run on a laptop if needed.

I've briefly tried to look around, but nobody seems to be asking this question with regards to simply parsing a BAM. Most memory-related topics for samtools seem to revolve around sorting.

So, is there a more resource-efficient way to parse BAMs progressively, or does the whole thing need to be decompressed into memory first (presumably that's what's happening) before the entries can be accessed sequentially?

Thanks!

Answer: It turns out that @i.sudbery was indeed correct. The memory-greedy bit was a dumb design choice in my script. I was storing information from each bam entry, only to count up the occurrences at the end. With 200000000 entries, the little bits of stored info added up...
My goal is much more efficiently achieved by starting all the counters right from the beginning and updating them entry by entry. This is really basic stuff, I'm not even sure why I didn't do it that way from the beginning.