cannin/biostars_qa huggingface dataset:
Question: Hi all.

I am using mpileup to extract some sequencing information. However, when I tried to run mpileup with below command

samtools mpileup -C 50 -vu -t DP -t DV -t DPR -t DP4 -f human_g1k_v37.fasta -l ms_TAR22_mpileup.bed TAR22_N_Filtered_Sorted_Makrdup_Readgroup_Realigned_Recal.bam

The results are nothing except head file of VCF.

I think that the problems occur due to the `-l` option because when I removed `-l` options, it works well but the file size was very large to handle. When I chose `-l` option and run it, as I expected, then it failed to generate sequencing information related to `-l` files. I don't know why mpileup command doesn't output anything but its head even though I followed the manual and ran with the latest samtools version(1.2).

For your reference, here is the data inside my bed files.

```
22 24829532 24829532 ADORA2A
18 12449779 12449779 SPIRE1
12 123005948 123005948 RSRC2
12 25398281 25398281 KRAS
12 22354743 22354743 ST8SIA1
7 93055822 93055822 CALCR
2 170374777 170374777 KLHL41
1 158532510 158532510 OR6P1
1 79002163 79002163 PTGFR
```

Help me to get it!

Answer: BED files are 0-based, half-open. What that means is that a SNV location needs to have a difference of 1 between the start and end location. In the case of your BED file, the locations appear to be in 1-based, closed coordinates. So, I suspect that mpileup is not using the BED file as you think. Try subtracting 1 from each element in column 2. Finally, make sure that the chromosome names match (including the "chr" or lack thereof).
Question: Hi all,

I wonder if there is any tool (in samtools ?) to quickly know whether a BAM file has been sorted or not.

(longer story: I've merged a set of sorted bams using 'samtools merge' and I've then used MarkDuplicates with the option: 'assume_sorted=true'. Am I right ?)

Thanks,

Pierre

UPDATE: I wrote a post about it : http://plindenbaum.blogspot.com/2011/02/testing-if-bam-file-is-sorted-using.html

Answer: You can use the sort order (SO) flag in the header to check if the file has been sorted:

```
% samtools view -H 5_110118_FC62VT6AAXX-hg18-unsort.bam
@HD VN:1.0 SO:unsorted
% samtools view -H 5_110118_FC62VT6AAXX-hg18-sort.bam
@HD VN:1.0 SO:coordinate
```

Unfortunately `samtools index` will work on both types of files without raising a command-line error, but will produce a smaller index file which is indicative of a problem:

```
% samtools index 5_110118_FC62VT6AAXX-hg18-sort.bam
% echo $?
0
% samtools index 5_110118_FC62VT6AAXX-hg18-unsort.bam
% echo $?
0
% ls -lh *.bai
chapman users 6.2M 2011-02-01 06:46 5_110118_FC62VT6AAXX-hg18-sort.bam.bai
chapman users 408 2011-02-01 06:47 5_110118_FC62VT6AAXX-hg18-unsort.bam.bai
```

**Edit:**

From the comments it appears as if samtools will sometimes raise an error in this case. Here is what an unsorted BAM file that finishes with a truncated index and no error messages looks like:

```
% samtools view 5_110118_FC62VT6AAXX-hg18-unsort.bam | head
HWI-EAS264:5:1:1353:936#0 77 * 0 0 * * 0 0 NTTTCCATTCCATTATGGATGATTCCATTCCATTGTATTC ######################################## XM:i:0
HWI-EAS264:5:1:1353:936#0 141 * 0 0 * * 0 0 TGGAACGGAATGGAATGGAATGGAATGGAATCAATCCGAC DDGGDEGGGGDG@BGGG=@G:EFFEEB4FF3BBB@BCB XM:i:0
HWI-EAS264:5:1:1691:931#0 99 chr18 28072372 255 40M = 28072771 439 NCTGGGTTTTTATGTCTTCCAGCTAGTAAGTGCTGGAGCT #322244442<:22<@@C@C@@;@;646666446664536 XA:i:1 MD:Z:0G39 NM:i:1
HWI-EAS264:5:1:1691:931#0 147 chr18 28072771 255 40M = 28072372 -439 GTTCTCCCGCCTGCACCTCTCTCTGTGTGTTGGCTTTGTT >B@D>ADDBDD>BD-DDDEDGGFHIIFI@HIIIHHIIIIH XA:i:0 MD:Z:40 NM:i:0
```
Question: I am new to bioinformatics in general and have the following situation:

I have just got 4 samples sequenced by Whole Exome Sequencing using Macrogen. I have received called variants from their bioinformatics pipeline as well as putting the generated fastq data from each sample through my own pipeline. I am looking to find a novel disease causing mutation in chromosome 2.

Amongst both pipelines a mutation a potentially interesting gene is called. However, the chromosomal coordinates are different. I understand that this is a somewhat probable situation but **is there a way that I might infer that one location is more likely than the other**?

I used GATK for variant calling and was thinking of using their Variant Quality Score Recalibration (VQSR) algorithms, but as this I believe, depends on previously determined SNPs I think is would bias against novel mutations.

Any ideas or suggestions would be appreciated.

Answer: You'll want to use the VQSR, since it'll decrease the false-positive rate. It's not so much that this biases against novel calls, but rather it uses the information gleaned from known sites to better gauge what's required for a real call.

For actual validation, you want to use an orthogonal technology (e.g., Sanger sequencing would suffice for a single gene).
Question: I am running `bedmap` and my output has many zeros after the decimal. I don't expect to have fractional counts. I looked at the bedmap help page and am just wanting clarification. Will the `--prec` flag allow me to remove the decimal and everything after it? I was considering using `--prec 0`.

Here is the code I am running:

./bedmap --fraction-map 0.5 --ec --echo --skip-unmapped --sum --delim '\t' hg38_chromSizes_500bpBins.bed $FILEIN > $OUT

And an example of the output:

chr1 10000 10500 28124.000000
chr1 10500 11000 2950.000000
chr1 11000 11500 0.000000
chr1 11500 12000 0.000000

Answer: Yes, that is what `--prec` is for. When processing signals, they are treated as floating point values internally, even if one starts with integers. If you expect output signal to be integers, it can be useful to specify `--prec 0`. See: https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#formatting-score-output

The `--skip-unmapped` option simply skips over any intervals in `hg38_chromSizes_500bpBins.bed` which do not overlap an interval in `$FILEIN`. If you expect that all intervals (bins) in `hg38_chromSizes_500bpBins.bed` are covered, then you do not need `--skip-unmapped`. However, it is often best to use this when using signal operations like `--sum`, etc. so that you skip intervals that would have no signal reported.

The `--echo` option will report intervals (bins) from `hg38_chromSizes_500bpBins.bed`, along with the results from `--sum`. If you want a BED file as output, then you would want to use `--echo` or one of the other interval-specific `--echo-*` options.
Question: Hi, I got stuck with an issue in analyzing my RNA-seq data. I have three tratments (Control, D1 and D2), and within each tratment I calculate DE genes between two groups (G1 and G2). Then I identify genes that become DE in both D1 and D2 compared to the control. Summarizing, I have (see also Figure):

- 15,000 genes representing the entire transcriptome
- 1,000 genes representing those that become DE as effect of treatment, filtering
those that are DE in the control. This group is made as union of D1
and D2 DE genes, respectively 450 and 600 genes.
- Intersecting the DE genes, I found that 50 of them are common to both
D1 and D2, 400 are specific to D1 and 550 are specific to D2.

![enter image description here][1]

Now I want to understand if this overlap is due by chance or not. So I'd like to have a feedback on my solution and, if wrong, a possible answer to my problem.

I have my observation: out of 1,000 DE genes common to both D1 and D2, 50 of them overlap. Then, I drawn 100 genes from both D1 and D2, and I counted how many overlap (bootstrapped 1,000 times). Let's say that the result of bootstrapped drawing is 10. Finally, I tested if the overlap observed and the one got by random drawing is different through a hypergeometric distribution.

I've done everything with R, using `phyper` as follows:

Here I have to use 4 values:

x = 10 (overlapping genes from drawing)
m = 50 (overlapping observed)
n = 950 (total DE genes - m)
k = 100 (number of genes drawn from D1 and D2)
lower.tail=FALSE

I got `0.007` which, if I understood correctly and my procedure is correct as well, means that the overlap I observe between D1 and D2 is not due by chance.

Any opinion about my approach?

Thanks!

[1]: https://i.stack.imgur.com/st2Kq.jpg

Answer: m = matrix(c(10, 550, 400, 15000-(10+550+400)), nrow=2)
fisher.test(m, alternative='greater')

You need to get the rates in each group from the total number of genes tested, so 15000 rather than 1000. You're wanting to know if the overlap is greater than expected by chance, thus the `alternative='greater'`. By chance you expect `450/15000 * 600 = 18` overlapping genes, which is more than you actually see. So since the p-value is near 1 you can't reject that you're seeing the overlap by chance (since you did a one-sided test, you can't then turn around and check the other one-sided test, since that'd be p-hacking).
Question: The following SNPs are from the GRCh38 build

1 158766655 158766655 C T
1 162023402 162023402 T C
1 169793666 169793666 T G
1 6554475 6554475 G A
3 179604597 179604597 A C
4 39062542 39062542 A C
5 177089630 177089630 G A
8 142786489 142786489 C T
8 142879785 142879785 C T
10 86654812 86654812 C T

how can I map them to GRCh37 (where possible)?

Thanks

Answer: Use NCBI's [remap][1] or UCSC's [LiftOver][2].

[1]: https://www.ncbi.nlm.nih.gov/genome/tools/remap
[2]: https://genome.ucsc.edu/cgi-bin/hgLiftOver
Question: Hi there!

I got a list of proteins ([UniProt ][1]ID) from a proteomics experiment and i am trying to build a model containing all the routes where my proteins are implicated. To achieve this, my first idea is to extract the [Reactome][2] route identifier from all the Uniprot IDs and then obtain the SBML for each one of those routes.

Joining the models (SBMLs) is a problem for future me, but i actually dont know how to obtain the SBML from [Reactome][2] with a list of identifiers like this (R-HSA-112409). I obtained the list of identifiers with a small python script.

I'm trying to build the model, not to get the routes more represented within my list of proteins. [Reactome][2] provides this type of analysis, but the SBML resulting from the top routes does not contain the same information that those accesed from the UniProt identifiers.

This is my first contact with SBMLs, so i would like to ask for advice for constructing the model (other databases providing SBMLs given a list of proteins for example), and if it were possible, for a way to obtain the SBMLs from [Reactome][2] using the route identifiers!

Thank you for your answers!

[1]: http://www.uniprot.org/
[2]: http://www.reactome.org/

Answer: Hi, when you start to work with SBML more and have questions, please feel free to post them on the [SBML forums][1] or [email the SBML Team directly][2]. The team has set up monitors on Biostars.org for SBML-related questions, but we don't tend to subscribe to Biostars directly simply because it seems a bit more bioinformatics-focused. Consequently, we might miss questions posted here or not find out about them until long after.

[1]: http://sbml.org/Forums
[2]: http://mailto:sbml-team@caltech.edu
Question: Hi all,

Does anybody know if there is a restriction for a read name length (say up to 10k characters)? Have you ever experienced that any software (like samtools, bwa, STAR, etc) complains about long read names?

Thanks

Answer: in htslib / sam.h https://github.com/samtools/htslib/blob/develop/htslib/sam.h

the read name is coded by

/* @field l_qname length of the query name */
uint8_t l_qname;

which is an unsigned byte. So the read name shouldn't be longer than 2^8=256.

it's even shorter in the read-sam method:

https://github.com/samtools/htslib/blob/f2150106bafda80ccd97970637091bb2799bc426/sam.c#L1305

_parse_err(p - q > 252, "query name too long");
Question: I did ggplot2 graph on R, the data is combined 2 dataframes, the plot has many 2 columns with same color each column representing each dataframe, the problem is each 2 columns for 2 dataframes has same colors but outlined with thin color, how to shade the data for one dataframe so that both same colors but one shaded one not shaded ?

```r
df$Cell_lines = 'cell_one'
data2$Cell_lines = 'cell_two'
df3=rbind(df,data2)
ggplot(df3, aes(x=Variant_Classification, y=Number_of_Genes,col =Cell_lines,fill = Variant_Classification)) +
geom_bar(stat="identity",position = 'dodge') + theme_minimal() + labs(title ="cell_one vs cell_two")+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
```

Answer: This code change the color and thickness of the border for better visualization:

df3 = rbind(df,data2)

ggplot(df3,aes(x=Variant_Classification,y=Number_of_Genes,col=Cell_lines,fill=Variant_Classification)) +
geom_bar(stat="identity",position='dodge',size =1)+
theme_minimal()+
labs(title="cell_twovscell_one")+
theme(axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank()) +
scale_color_manual(values = c("cell_one" = "red","cell_two"="black"))
Question: Hello,

I have about 30 VCF files with file names as ID_001.new.vcf. I want to extract only the "ID_001" part from the file name and change it in the header line of the VCF file where "Sample1" is given.

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1

So that the result looks like that:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ID_001

How can I do it ? I tried to use echo in bash and extract the IDs from the Filename but I am unable to iterate it to change inside the file.
Thanks for your help.

Answer: In bash this should do.

for i in *.new.vcf
do
ID_NAME=$(basename "$i" .new.vcf)
sed -i "1s|Sample1|$ID_NAME|g" $i
done

Caution: I have used -i with sed. So the actual files will get edited in place.

Now added 1s also as to limit the replacement to first line alone.
Question: A philosophical question, but one that has been posed to me by two independent PIs in the last month. If there are two biological replicates for the same ChIP sample of a ChIP-seq experiment and one of them performs "better" than the other (in terms of number and quality of detected regions) should they be combined?

My own personal take on this is that they should not be combined. Peak calling should be performed as replicates and then the intersect between the replicates taken. But I would be interested to see if anyone can see merit in combining biological replicates.

Thanks.

Answer: The first thing I do prior to any of this is, post-alignment, run some QC using the ENCODE tools, in particular the [phantompeak tools][1] strand cross-correlation (NSC, RSC, and to some degree QTag) analysis and PCR bottleneck metrics. If you have biological reps then [IDR][2] is a very good idea as it checks consistency between replicates. [CHANCE][3] is also a highly recommended tool for QC. **Note**: IDR only works for TF and is not really developed for broad peaks (histone modification), so keep that in mind.

Do these on **all** your samples, including inputs; you should see no cross-correlation in inputs and decent cross-correlation on your IP. The level of the latter being dependent on whether it is narrow vs. broad-peak; if you are doing TF then use the criteria outlines in the phantompeak tools page. We have seen some cases with inputs where you will see some cross-correlation, possibly indicative of biases in fragmentation (open chromatin fragments more easily than closed, for instance), and these will influence your peak calls. If you have access to a Covaris I highly recommend it as we have found the best overall results from samples prepped from those, but they're quite pricy.

I also recommend if you are performing ChIP-Seq runs that ENCODE has previously performed (same organism, same factor) you can run a comparison against theirs to get an idea of how your samples come out (there is an on-line speadsheet [on the ENCODE site][4] with this information).

**Very useful**: Use IGV and calculate coverage metrics (TDF). Run some preliminary peak calls on your best sample and visually check peaks. Use ENCODE data sets from IGV to compare to if possible!

After all that, if the samples are quite comparable based on IDR and peak metrics then combining them is probably fine for peak-calling, but if they are different then you are bound by your worst sample and thus would likely be introducing noise into the peak calling (e.g. false positives).

[1]: https://code.google.com/p/phantompeakqualtools/
[2]: https://sites.google.com/site/anshulkundaje/projects/idr
[3]: http://songlab.ucsf.edu/CHANCE.html
[4]: https://genome.ucsc.edu/ENCODE/qualityMetrics.html
Question: I'm trying to use WGCNA for weighted network analysis. I'm getting trouble how to decide this parameter. It's clear how to define and why this has been defined.

#Choosing a soft-threshold to fit a scale-free topology to the network

powers = c(c(1:10), seq(from = 12, to=20, by=2));

What is the meaning of this powers and how to decide this value?

Answer: This parameter must be >1 and is simply the power to which the correlations are raised. This increases the contrast between high and low values, in effect a form of soft thresholding. In the [paper][1], the authors suggest setting it so that the resulting network has a scale-free topology. See the pickSoftThreshold() function in the [WGCNA R package][2].

[1]: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-9-559
[2]: https://CRAN.R-project.org/package=WGCNA
Question: Hi all,
Currently, after mapping, I got a sam file which contained approximately 30 million mapped reads. Then, by using bedtool bamtobed, I was able to obtain a bed file. However, there were so many many similar reads.
Therefore, in order to reduce the size of the final bed file, is there any way that I can mention each read only once as well as its number of presentation in the sam file.

----------

For instance:

----------

var_1 0 15 ATGCATGCATGCCGTA

----------
var_1 0 15 ATGCATGCATGCCGTA

----------
var_1 0 15 ATGCATGCATGCCGTA

----------
var_2 5 20 ATGCATGCGGGCCCC

----------
Will become:

----------
var_1 0 15 ATGCATGCATGCCGTA 3

----------

var_2 5 20 ATGCATGCGGGCCCC 1

Thank you in advance!

Answer: Use `uniq -c` your bed file.

uniq -c mybed.bed > uniqbed.bed

or

echo "var_1 0 38 ATGCATGCATGCCGTA
var_1 0 38 ATGCATGCATGCCGTA
var_1 0 38 ATGCATGCATGCCGTA
var_2 5 40 ATGCATGCGGGCCCC" | uniq -c

3 var_1 0 38 ATGCATGCATGCCGTA
1 var_2 5 40 ATGCATGCGGGCCCC
Question: I need to download WGS for healthy women, ideally sequenced on Illumina HiSeq 2000. I intitially looked into 1000G data ([http://www.internationalgenome.org/data/][1]), but it doesn't seem to specify which machine was used for sequencing. I am interested in low complexity sequences and thus worried the machine type could make a difference to my downstream analysis.

Thank you.

[1]: http://www.internationalgenome.org/data/

Answer: How about some of the [CEPH pedigree 1463 data][1] (Hiseq2000, 50x). For these, you should also find plenty of validated variant calls if you search around.

[1]: https://www.illumina.com/platinumgenomes.html
Question: Guys, I've wrote a script to extract sequences between the 23S rRNA and 16S rRNA gene in Python using Biopython. I can't figure out why, even if I've set explicit conditions, it keeps extracting randome sequences into the genome. Here's the code:

for f in glob.glob(cwd+'/genomes/*.gbk'):
gbank = SeqIO.parse(open(f,'r'), "genbank")
with open(output_file, "w") as nfh:
for genome in gbank:
for gene in genome.features:
if gene.type=='rRNA':
if gene.location.strand == -1: # wanted just the reverse strand
if 'product' in gene.qualifiers:
if '16S' in gene.qualifiers['product'][0]:
end = gene.location.nofuzzy_start
elif '23S' in gene.qualifiers['product'][0]:
start = gene.location.nofuzzy_end
if abs(start - end) < 50000:
intergenic_sequence.append(">%s|ITS_full_sequence|%s\n%s" % (genome.description, gene.location.strand, genome.seq[start:end]))

I've set `gene.type` to `rRNA` and it extract a sequence into a random CDS, nothing to do with ribosomal RNA. This problem happens with several genomes tested (RefSeq records, well annotated).

Answer: Many many thanks for this solution goes to @jrj.healey . Apparently the logic of features selection is the same, so I suddenly don't realize why this selection works instead of the previous.

for f in glob.glob(cwd+'/genomes/*.gbk'):
gbank = SeqIO.parse(open(f,'r'), "genbank")
with open(output_file, "w") as nfh:
for genome in gbank:
# trait for ITS located onto reverse strand
rnas_rev = [feat for feat in genome.features if feat.type == 'rRNA' and feat.location.strand == -1]
s_rev = [ r for r in rnas_rev if '16S' in r.qualifiers['product'][0]]
lo_rev = [ r for r in rnas_rev if '23S' in r.qualifiers['product'][0]]
# trait for ITS located onto forward strand
rnas_fw = [ feat for feat in genome.features if feat.type == 'rRNA' and feat.location.strand == -1]
s_fw = [ r for r in rnas_fw if '16S' in r.qualifiers['product'][0]]
lo_fw = [ r for r in rnas_fw if '23S' in r.qualifiers['product'][0]]
# i select 16S genes
# j select 23S
for i,j in zip(s_rev,lo_rev):
....
for i,j in zip(s_fw,lo_fw):
....
Question: Hi everyone,

I have a doubt regarding how DESeq2 handles the HTSeq counts assigned to special classes when performing normalization & differential gene expression.

Apart from actual genes, HTSeq assigns reads to the five classes below. I have an approximation of the mean of the counts assigned to each class for my dataset:

1. `no_feature` (~40112886)
2. `ambiguous` (~9732)
3. `too_low_aQual` (0)
4. `not_aligned` (0)
5. `alignment_not_unique` (~4294028)

Now, I am going to use DESeq2 to normalize and get differentially expressed genes for my data. So how does DESeq2 handle these classes? Does it remove them during normalization, uses them in the normalization process or do we have to remove these rows manually before the normalization?

Thanks!

Answer: It removes them when making the DeseqDataSet object that's used in the subsequent steps.

**Edit**: If you're curious, here are the relevant lines of `DESeqDataSetFromHTSeqCount`:

oldSpecialNames <- c( "no_feature", "ambiguous",
"too_low_aQual", "not_aligned",
"alignment_not_unique" )
# either starts with two underscores
# or is one of the old special names (htseq-count backward compatability)
specialRows <- (substr(rownames(tbl),1,2) == "__") |
rownames(tbl) %in% oldSpecialNames
tbl <- tbl[ !specialRows, ]

**Edit2**: I should mention that if you use one of the other functions, such as `DESeqDataSetFromMatrix` or `DESeqDataSet`, you need remove these fields yourself.
Question: Hello every one,

It's about a month that I stuck with picard. Every command that I try to run with picard failed.
And for all of them I get the same error:

> ERROR: Invalid argument '-I' , '-O'

When I try this argument I= or R=

I get this:

> NOTE: Picard's command line syntax is changing.
>
> ** For more information, please see: **
> https://github.com/broadinstitute/picard/wiki/Command-Line-Syntax-Transition-For-Users-(Pre-Transition)
>
> ** The command line looks like this in the new syntax:
> ** MergeVcfs -I sample1.vcf -I sample2.vcf.vcf -O Output.vcf.gz

can anyone tell me if I am doing sth wrong or it's a bug in picard?

I really appreciate any help!

Answer: works on my machine using https://github.com/broadinstitute/picard/releases/download/2.18.20/picard.jar

java -jar /path/to/picard.jar MergeVcfs I=/path/to/S1.vcf.gz I=/path/to/S1.vcf.gz O=/path/to/out.vcf.gz
(...)
[Sun Dec 16 14:57:17 CET 2018] picard.vcf.MergeVcfs done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=40632320

yes there is the message ` NOTE: Picard's command line syntax is changing.` but that's just an information.