cannin/biostars_qa huggingface dataset:
Question: Hello BioStars,

Is there a tool which generates mapping quality score distributions from bam files?

I know I could potentially do this myself, but I am looking for something which would essentially do the work for me, and also something which plugs in nicely to multiqc

I have had a quick search, and it is surprisingly difficult to find a tool which provides this functionality, so I am just wondering if I am missing something obvious here? The obvious candidate would be something like picard, but picard does not seem to offer this functionality from what I can tell

Thanks in advance!

Answer: I'm pretty sure that [Qualimap BamQC][1] (not to be confused with the BamQC tool from Babraham mentioned above) calculates and plots this. Qualimap is supported by MultiQC, however off the top of my head, I don't think that MultiQC replicates this plot in reports. I doubt it would be particularly difficult to add if you fancy it (or opening a GitHub issue with some example data), or you could probably sidestep that by using MultiQC's "[custom content][2]" feature.

[1]: http://qualimap.conesalab.org/doc_html/analysis.html
[2]: https://multiqc.info/docs/#custom-content
Question: Hi all,

I have fasta sequence of some proteins and I want to convert fasta format to phylip file format to build phylogenetic tree using ggtree. I tried online EMBOSS seqret tool to convert fasta file to phylip format but I got error when i read in ggtree.

my input sequence

>proteinsA
MGDSRDLCPHLDSIGEVTKEDLLLKSKGTCQSCGVTGPNLWACLQVACPYVGCGESFADH
RTDKKPALCKSYQKLVSEVWHKKRPSYVVP
>proteinsB
MTGSNSHITILTLKVLPHFESLGKQEKIPNKMSAFRNHCPHLDSVGEITKEDLIQKSLGT
SHVSFP
>proteinsC
MGDSRDLCPHLDSIGEVTKEDLLLKSKGTCQSCGVTGPNLWACLQVACPYVGCGESFADD
ITTEETMEEDKSQSDVDFQSCESCSNSDRAENENGSRCFSEDNNETTMLIQDDENN

and EMBOSS seqret output is..

3 116

proteinsA MGDSRDLCPH LDSIGEVTKE DLLLKSKGTC QSCGVTGPNL WACLQVACPY
proteinsB MTGSNSHITI LTLKVLPHFE SLGKQEKIPN KMSAFRNHCP HLDSVGEITK
proteinsC MGDSRDLCPH LDSIGEVTKE DLLLKSKGTC QSCGVTGPNL WACLQVACPY

VGCGESFADH RTDKKPALCK SYQKLVSEVW HKKRPSYVVP ----------
EDLIQKSLGT SHVSFP---- ---------- ---------- ----------
VGCGESFADD ITTEETMEED KSQSDVDFQS CESCSNSDRA ENENGSRCFS

---------- ------
---------- ------
EDNNETTMLI QDDENN

But I got error in reading this phylip file...

tree <- read.phylip("emboss_seqret_output.txt")

Error in read.phylip("emboss_seqret_output.txt") :
input file is not phylip tree format...

Could you please help me what is problem with my input file or can you please suggest me to some alternative ways.

Thanks a lot.

Answer: I found a nice tutorial to build Phylogenetic Trees from fasta sequence...

http://www.cbs.dtu.dk/courses/biosys/binfintro/phylogeny.php

Step 1: Open the sequence file (fasta), select the entire file, and copy the sequences.

Step 2 : Align the sequences in using the mafft server at EBI with default settings as follows

Step 3: Open the TreeHugger web server. (The TreeHugger server constructs a neighbor joining tree from an aligned set of sequences).

Step 4: Download data in Newick/Phylip format" (treehugger_newick.nwk)

Step 5: Visualizing using ggtree

library(ggtree)
tree <- read.tree("treehugger_newick.nwk")
ggtree(tree) + geom_tiplab()
Question: Hi, I'm wondering what exactly is the meaning of an aligner being "splice aware". I know it has to do with the mapping of reads spanning splice junctions, but as someone pretty new to RNA-seq and molecular biology, that's not quite enough for me to grasp the concept.

The following is my best reasoning of its meaning. RNA-seq reads are derived from mature mRNA, so there's typically no introns in the sequence. But aligners use a reference genome to aid in the process, so a read spans (what in the actual transcript are) two exons, while the reference would have one exon followed by an intron. So the reference genome would find a matching sequence in only one of the exons, while the rest of the read would not match the intron in the reference, so the read can't be properly aligned. A splice-aware aligner would know not to try to align RNA-seq reads to introns, and would somehow identify possible downstream exons and try to align to those instead, ignoring introns altogether.

Is this anywhere close to the meaning of splice-aware? And if so, would a splice-unaware aligner properly align RNA-seq data, given a reference transcriptome?

Answer: I think, that pretty much hits the point. :) The major problem is that introns not only vary in length but that they can also be very long. An DNA-DNA aligner (i.e., splice-unaware) would have to introduce a long gap in the mapping of a read to span an intron. This is not desired for DNA read mapping and might lead to false mappings. Without transcript information/restriction it is quite likely that the aligner is able to find a genome-sequence that matches the remaining read sequence. But that can be at any position downstream...

Regarding your second question, the splice-aware Tophat internally uses the splice-**unaware** Bowtie. So yes, in principle it is possible to use a splice-unaware aligner for mapping RNA-Seq reads to a transcriptome. However, using dedicated *RNA-Seq read mapper* or *RNA-Seq read mapping work flow* might give you better results by taking care of possible caveats etc. (read the papers for this ;-) ).

PS: See also https://www.biostars.org/p/175068/
Question: Hello,
I am having a hard time opening a very large bedgraph. I have been suggested to split my bam file by chromosome with chrom-bed.py but it didn't work.
Is there any other alternative?
Thanks,
GP.

Answer: Try samtools: `samtools view -?`

A region should be presented in one of the following formats:
`chr1',`chr2:1,000' and `chr3:1000-2,000'. When a region is
specified, the input alignment file must be an indexed BAM file.

something like `samtools view in.bam chr1 > chr1.bam` should work
Question: I have some genomic coordinates and I want to retrieve their RNA sequence especially the ones for reverse complement ones (since RNA sequence would not be strand specific). Where could I find such tool or package?

Answer: Hello,

you can also go to [the UCSC Table Browser][1], upload your input file as a custom track ("add custom track" button), select your genome, choose "sequence" as output format, press "get output" and choose "mRNA" from the following menu. This would output your mRNA sequences in FASTA format.

[1]: https://genome.ucsc.edu/cgi-bin/hgTables
Question: Is there a way to use custom matrix for nucleotide aligner in Biopython pairwise2 module?
I want to decrease penalties for A->G, C->T substitutions.

import Bio
from Bio import pairwise2

seq1='AGTCGATCATCGGA'
seq2='ATTAGATCATCGGGGA'

alignment = pairwise2.align.localms(seq1, seq2, 2,-1,-2,-1,one_alignment_only=True)
print(format_alignment(*alignment[0]))

AGTCGATCATCGGA--
|.|.|||||||||
ATTAGATCATCGGGGA
Score=20

Answer: how about trying somethig like this

import Bio
from Bio import pairwise2
from Bio.pairwise2 import format_alignment

seq1 = 'AGTCGATCATCGGA'
seq2 = 'ATTAGATCATCGGGGA'

match_dic = {}
match_dic[('A', 'A')] = 2
match_dic[('A', 'C')] = -1
match_dic[('A', 'T')] = -1
match_dic[('A', 'G')] = 1
match_dic[('C', 'A')] = -1
match_dic[('C', 'C')] = 2
match_dic[('C', 'T')] = 1
match_dic[('C', 'G')] = -1
match_dic[('G', 'A')] = -1
match_dic[('G', 'C')] = -1
match_dic[('G', 'T')] = -1
match_dic[('G', 'G')] = 2
match_dic[('T', 'A')] = -1
match_dic[('T', 'C')] = -1
match_dic[('T', 'T')] = 2
match_dic[('T', 'G')] = -1

alignment = pairwise2.align.localds(
seq1, seq2, match_dic, -2, -1, one_alignment_only=True)
print(alignment[0])
print(format_alignment(*alignment[0]))
Question: Hey guys,

I'm having issues trying to run a multivariate linear mixed model (lmm) Genome Wide Association Study (GWAS).

My input files & genotyping files are fine, but when I tried to run the GWAS in [Gemma][1], the program stalled.

I've no issue running GWAS on individual phenotypes, but I've been told that a multivariate analysis using Gemma won't work; because my input phenotype values are biomodal & apparently Gemma can only run multivariate GWAS on data with a normal distribution.

Here is the input I used on Gemma:

gemma -bfile multivariate_genotypes -k output/kinship.cXX.txt -n 2 3 -lmm -o multivate_PC1+PC2_lmm_20JUN2017

Please let me know if I'm using Gemma incorrectly or if you can suggest more appropriate software.

Thanks,

Ronan

[1]: http://www.xzlab.org/software.html

Answer: It is not as much GEMMA but possibly the linear model itself. GEMMA's manual (4.2.3?) has a reference to the linear model's robustness to model misspecification but also to the ability to model binary phenotypes and the authors provide a further explanation in

Xiang Zhou, Peter Carbonetto, and Matthew Stephens. Polygenic modelling with Bayesian sparse linear mixed models. *PLoS Genetics,* 9:e1003264, 2013.

I am not sure how robust the multivariate version would be to model misspecification - even if you analyzed your data in a way that would represent the bimodal distribution as a categorical (using data-driven approach - maybe even latent cluster/profile modeling - that would also answer questions about the multivariate dimensionality of your data?) /dichotomous DV value. However, that is, partly, an empirical question - I would run both analyses using the linear model. Why is the mixed aspect of it of particular interest as opposed, to, say, using categorical regression with PC's as covariates? Due to relatedness?

Anyway, do share your results! (I realize this question is 11 months old).
Question: I can extract reads mapped to a single reference sequence by first aligning the reads using, say, BWA, and processing with SAMtools:

samtools view -b -f 0x2 alignment.sam | samtools fastq - -1 mapped_1.fastq -2 mapped_2.fastq

Is there a similar pipeline I can use to extract reads mapped to multiple references? That is, if a read was found in any of the given references, it should be included in the output.

Here's an idea I have:

1. align the reads to each reference separately (producing
multiple SAM files)
2. extract the reads mapped to each
reference using the above pipeline
3. `sort` and `uniq` the names of the reads from all output files, storing the names in file.txt
4. `grep` the names in file.txt from one of the SAM files and store in a new SAM file
5. convert SAM to FASTQ

I'm guessing this is probably not the most efficient way to solve the problem. It also doesn't catch when a read might be mapped to one of the references without its pair. In this case, I'd like to extract both reads.

Answer: The most efficient way would be to map to all the references at once. You can do this with any aligner if you just concatenate the references. Alternatively you could map to all at once using BBSplit, which can produce one fastq output file per reference, containing all the reads (or pairs) that best map to that reference. But if you simply want all reads that map to any reference in a single fastq file, along with their unmapped mates, you can do this:

1) Concatenate the references with cat.
2) Map the reads with this command:
bbmap.sh ref=references.fa in1=read1.fq.gz in2=read2.fq.gz outm=mapped.fq.gz outu=unmapped.fq.gz

If one or both reads in a pair map to anything, both will go to "outm"; if neither map, both will go to "outu". The output files will be interleaved, but if you use "outm1" and "outm2" (and etc.) you can get the output in dual files.
Question: I would like to reshape a dataframe and repeat some informations about Name strings ("Pro1", "Pro2" and "Pro3").

My input:

Name Chr Position NE001 NE002 NE003
Pro1 1 25 -0.3 0.2 0.4
Pro2 2 23 -0.4 0.2 -0.3
Pro3 3 30 -0.3 0.2 0.4

My expected output:

Sample Pro1 Chr Position Pro2 Chr Position Pro3 Chr Position
NE001 -0.3 1 25 -0.4 2 23 -0.3 3 30
NE002 0.2 1 25 0.2 2 23 0.2 3 30
NE003 0.4 1 25 -0.3 2 23 0.4 3 30

Cheers!

Answer: Given your input as a dataframe or matrix named `d`:

d2 <- t(d[,-c(1:3)])
d2 <- cbind(d2, matrix(rep(d[,2], each=nrow(d2)), nrow=nrow(d2)),
matrix(rep(d[,3], each=nrow(d2)), nrow=nrow(d2)))
colnames(d2) <- c(as.character(d$Name),rep("Chr", nrow(d2)), rep("Position", nrow(d2)))
d2 <- d2[,c(1,4,7,2,5,8,3,6,9)]

The output is:

> d2
Pro1 Chr Position Pro2 Chr Position Pro3 Chr Position
NE001 -0.3 1 25 -0.4 2 23 -0.5 3 30
NE002 0.2 1 25 0.2 2 23 0.2 3 30
NE003 0.4 1 25 -0.3 2 23 0.4 3 30
Question: I have WGS paired-end data and I want to use BreakDancer to call structural variants.

Both my case and controls bam files do not have, in their header, the LB field in @RG tag and BreakDancer reports this lack. I can re-header bam files and add the missing LB easily, however I cannot understand how this information is important for BreakDancer and which are the recommended criteria to add them, i.e. a unique LB for all samples or specific for each of them? Is it referred to the Library Preparation Kit?

Could you provide me some hints to solve the problem? Thank you

Answer: Just use picard's AddOrReplaceReadGroups tool to replace the read groups. Assuming that each sample was sequenced only once, just set the library (LB) field to be the same as the sample (SM). So java -jar AddOrReplaceReadGroups.jar I=original.bam O=new.bam SM=sample1 LB=sample1 ....
Question: I can't find a compresensive guide to do GO term enrichment. Is there an accepted tutorial or review paper?

Thanks!

Answer: http://www.bigre.ulb.ac.be/courses/statistics_bioinformatics/practicals/ASG1_2012/go_statistics_td/go_statistics_td.html

http://mengnote.blogspot.com/2012/12/calculate-correct-hypergeometric-p.html

http://users.unimi.it/marray/2007/material/day4/Lecture7.pdf
Question: Hi,

I just had a quick question that I started to confuse myself with. When I am about to index my reference sequence for some RNA-Seq data and I got confused if I use the cDNA or CDS or transcript fasta file.

In the past I used a transcript fasta but I can't seem to find one for the species an am working with from ensembl (http://useast.ensembl.org/info/data/ftp/index.html/)

Maybe I can just use the whole genome for the reference? I am using Salmon to map the reads.

Thanks in advance

Answer: You should not map against the genome using Salmon. You can either download a transcriptome file, or a genome file and transcript annotations, and use a tool like gffread to extract the transcript sequences. You most likely want to quantify against the cDNA to account for features such as UTRs.
Question: Dear fellows:

I have some ATAC bam files(A549 in drugged and control condition) which have different sequence depths that need to convert to bigwig.

The tools I use for conversion is **deeptools bamCoverage**, and I use **deeptools multiBamSummary** to generate scale factor for each bam file. then I use **bamCoverage** with `--normalizeUsing RPGC` and `--scaleFactor of each bam` file to generate bigwig.

After I got the bigwig file of each sample, I load them into IGV. My propose is to display different heights of peak pileups(eg. higher peak pileups in treated samples and lower pileups on control) at target differential access positions(already known) between those samples.

Here are my questions:
1. Does my normalization way(use `RPGC` and `scale factors` when converting bam to bigwig by **bamCoverage**) make sense?
2. when I load **already normalized** bigwig files into IGV without adjusting any IGV params, does the displayed track height of each sample could give the expected trend(the higher track means actual reads pileups and vice versa), or do I still need to use IGV's function **autoscale** or adjust **data range** to some range to let the peak fit the trend?

Here are two figures related to my question: Fig. A and Fig.B both use normalized bigwig tracks by `RPGC` and a `scaling factor` of positive and negative treated and control samples. the red rectangle area means the expected differential access positions will show the different track heights.

In Fig A, I use IGV autoscale function to 3 tracks; In Fig B, I set the data range of 3 tracks both to 50. The position 1 in two figures gives the conflict interpretation: in Fig A auto-scaled track, their height seems the same; but in Fig B tracks be set the same data range, the negative_treated track has the highest peak.

##### Fig. A: IGV auto scaled
![ABCC3_autoscale][1]

##### Fig. B: IGV set the same data range to 50
![ABCC3_same_range][2]

Can I have your advice on solving my issues?

Many thanks!

[1]: /media/images/3fd7d10e-a8e8-48b5-ab29-474b2f94
[2]: /media/images/9957d964-330b-4387-aa7a-230a2e51

Answer: Yes, the RPGC normalization make sense.

Considering the visualization within IGV, I would recommend to using either option #2 or the 'group autoscale' feature (select all tracks, right click, left click on `group autoscale`). With that setting, all tracks will have the same range, like in your second example except that the range will automatically fit all the data. In both option #2 and , higher peak means more coverage which is the correct interpretation.
Question: Hello,

After searching in this user group and analyzing the output of picard MarkDuplicates, it seems that picard does not output which of the reads it determined to be optical duplicates. That is, picard gives the number of optical and total duplicates in their metric file, but it only marks in its output bam file which reads are duplicates (optical or otherwise), and not which of those are optical duplicates.

I was requested to find a way to remove the optical duplicates and the only way I can think of doing this is to change the picard source code itself to output such a flag into the bam file (which could then be used to easily filter on). Before embarking on this task however, I'd like to know if someone has already done this and would this person be so kind as to share the code to spare me the trouble?

Thanks a lot!!!

Anna

Answer: I contacted the Broad Institute, and Nils Homer (samtools team) informed me that they already knew about one bug relating to this. I tested the fix, and this bug affected not only the number of optical duplicates but also the estimated library size statistic. To get the bug fix you need version >= Picard 1.120.

Unfortunately I still think there are additional bug(s), and I'm working with Nils on this. I'll let you know what else I find out!!
Question: Hi all,

I am using clusterProfiler to perform KEGG enrichment. I have a list of gene Symbol (Rat genes). First I need to translate the symbols to EntrezID.

I installed the Genome wide annotation for Rat:

source("https://bioconductor.org/biocLite.R")
biocLite("org.Rn.eg.db")

Then I tried to run the command:

x <- c("GPX3", "GLRX", "LBP", "CRYAB", "DEFB1", "HCLS1", "SOD2", "HSPA2",
"ORM1", "IGFBP1", "PTHLH", "GPC3", "IGFBP3","TOB1", "MITF", "NDRG1",
"NR1H4", "FGFR3", "PVR", "IL6", "PTPRM", "ERBB2", "NID2", "LAMB1",
"COMP", "PLS3", "MCAM", "SPP1", "LAMC1", "COL4A2", "COL4A1", "MYOC",
"ANXA4", "TFPI2", "CST6", "SLPI", "TIMP2", "CPM", "GGT1", "NNMT",
"MAL", "EEF1A2", "HGD", "TCN2", "CDA", "PCCA", "CRYM", "PDXK",
"STC1", "WARS", "HMOX1", "FXYD2", "RBP4", "SLC6A12", "KDELR3", "ITM2B")

eg = bitr(x, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Rn.eg.db")

Then I got the warning message:

> 'select()' returned 1:1 mapping between keys and columns Warning
> message: In bitr(x4, fromType = "SYMBOL", toType = "ENTREZID", OrgDb =
> Rat) :
> 98.21% of input gene IDs are fail to map...

Could anyone tell me what happened and how should I do?

Many thanks,
Stanley

Answer: ```
> paste0(substring(x, 1, 1), tolower(substring(x, 2))) -> x2
> x2[x2=="Pvr"] = "PVR"
> eg = bitr(x2, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Rn.eg.db")
'select()' returned 1:1 mapping between keys and columns
> eg
SYMBOL ENTREZID
1 Gpx3 64317
2 Glrx 64045
3 Lbp 29469
4 Cryab 25420
5 Defb1 83687
6 Hcls1 288077
7 Sod2 24787
8 Hspa2 60460
9 Orm1 24614
10 Igfbp1 25685
11 Pthlh 24695
12 Gpc3 25236
13 Igfbp3 24484
14 Tob1 170842
15 Mitf 25094
16 Ndrg1 299923
17 Nr1h4 60351
18 Fgfr3 84489
19 PVR 25066
20 Il6 24498
21 Ptprm 29616
22 Erbb2 24337
23 Nid2 302248
24 Lamb1 298941
25 Comp 25304
26 Pls3 81748
27 Mcam 78967
28 Spp1 25353
29 Lamc1 117036
30 Col4a2 306628
31 Col4a1 290905
32 Myoc 81523
33 Anxa4 79124
34 Tfpi2 286926
35 Cst6 171096
36 Slpi 84386
37 Timp2 29543
38 Cpm 314855
39 Ggt1 116568
40 Nnmt 300691
41 Mal 25263
42 Eef1a2 24799
43 Hgd 360719
44 Tcn2 64365
45 Cda 362638
46 Pcca 687008
47 Crym 117024
48 Pdxk 83578
49 Stc1 81801
50 Wars 314442
51 Hmox1 24451
52 Fxyd2 29639
53 Rbp4 25703
54 Slc6a12 50676
55 Kdelr3 315131
56 Itm2b 290364
```

I think we should raise this issue to Bioconductor. IMO, Mapping ID with case insensitive is necessity.
Question: Hello everybody,
I need to obtain as much more as possible protein homology structures (.pdb) from a few hundred genes of Rattus norvegicus. Do someone know if is there a way to obtain this structures in an high-throughput manner?
Any local/web-based program or Python/R module would be great.

Answer: There is a good chance that models for your proteins may already exist:

- https://modbase.compbio.ucsf.edu/
- https://swissmodel.expasy.org/repository

Or you can use the links below to model them:

- https://modbase.compbio.ucsf.edu/modweb/
- https://swissmodel.expasy.org/
Question: Is there any resource of standardized and uniformly stylized (and potentially, but not necessarily, even free) icons of popular model organisms (including E. coli and other popular unicellular organisms)?

I wish to use it for making a small tabular-like panel with ~10 organisms, next to the main figure panel, which would show the detailed data for one organism.

Some model organisms are not yet included in unicode.

Answer: See the [PhyloPic][1] project.

Edit: and the [rphylopic][2] project to access PhyloPic from R.

[1]: http://phylopic.org
[2]: https://github.com/sckott/rphylopic
Question: Hello, I have a set of .vcf files with phased genotypes that I would like to convert to standard haplotype format (.hap). Is there any tool that makes this?

Here is a toy example of a .hap file of a population of 2 human individuals, with 12 SNPs each:

1 1 2 2 1 2 2 2 1 1 1 1 1
2 1 1 1 2 1 1 2 2 1 2 1 2
3 1 1 2 1 2 2 1 1 2 2 1 2
4 2 2 2 2 2 2 1 1 1 1 1 1

Each line represents a chromosome, the first element being the haplotype ID. 1 and 2 represent ancestral and derived alleles, respectively.

Many thanks

Answer: Assuming all of your genotypes are phased, this should be doable with recent plink 2.0 alpha builds (which you can download from https://www.cog-genomics.org/plink/2.0/ ):

plink2 --vcf phased.vcf --export haps --out new_filename_prefix
Question: FASTQC, for example, doesn't seem to have a publication associated with it. How would you cite it?

Answer: I suggest you use the author, date and link.