cannin/biostars_qa huggingface dataset:
Question: (I also asked this question on the Bioconductor support site [here][1].)

I'm stumped. I'm trying to plot a few transcripts at the same time, given transcript names and a TxDb. These are examples of approaches I've tried:

# ------------------------------------------------------------------------------
# Setup:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(Gviz)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

# ------------------------------------------------------------------------------
# Try 1:
gr <- GenomicFeatures::exons(
txdb,
vals = list(tx_name = c("uc001aaa.3", "uc010nxq.1")),
columns = list("EXONNAME", "TXNAME", "GENEID"))
track <- Gviz::GeneRegionTrack(gr)
Gviz::plotTracks(track)
# Creates a plot, but doesn't show transcript grouping

![enter image description here][2]

# ------------------------------------------------------------------------------
# Try 2
gr <- GenomicFeatures::transcripts(
txdb,
vals = list(tx_name = c("uc001aaa.3", "uc010nxq.1")),
columns = list("EXONNAME", "TXNAME", "GENEID"))
track <- Gviz::GeneRegionTrack(gr)
Gviz::plotTracks(track)
# Creates a plot, but has no exon/intron information

![enter image description here][3]

# ------------------------------------------------------------------------------
# Try 3
gr <- exonsBy(txdb, by = "tx", use.names=TRUE)[c("uc001aaa.3", "uc010nxq.1")]
track <- Gviz::GeneRegionTrack(gr)
# Error in .fillWithDefaults(DataFrame(chromosome = as.character(seqnames(range)), :
# Number of elements in argument 'feature' is invalid

None of these work for me. I want to display exon/intron-structures and have arrows between exons. Does anyone have suggestions?

[1]: https://support.bioconductor.org/p/80221/
[2]: http://s8.postimg.org/dv9mopwhh/try1.png
[3]: http://s23.postimg.org/hi1ta61pn/try2.png

Answer: I found a solution that worked! Typical to find it shortly after asking the question, but here it is:

# ------------------------------------------------------------------------------
# Try 4
gr <- exonsBy(txdb, by = "tx", use.names=TRUE)[c("uc001aaa.3", "uc010nxq.1")]
gr <- unlist(gr)
elementMetadata(gr)$transcript <- names(gr)
track <- Gviz::GeneRegionTrack(gr)
Gviz::plotTracks(track)

![enter image description here][1]

Using exonsBy gives a GrangesList with one GRanges object for each transcript. Unlisting the GrangesList creates a named (important!) GRanges object with all exons. Using those names to set the $transcript metadata column solved the problem.

[1]: http://s16.postimg.org/ong5l77w5/try4.png
Question: I have a dataframe in R and I want to plot a subset of the plot as a line graph in ggplot. I have 8 different variables, with no guarantee all 8 will appear in the subset I want to plot. How can I for ggplot to assign variable A to a particular color code #B35806 and H to #542788? I tried to assign this to the dataframe itself (a column where if A is present, #B35806 would be) and calling on that in ggplot but that did not help

geom_line(aes(color=Line_color))

Advice would be very useful to this R noob.

Answer: I found a solution

p + scale_colour_manual(values = c("A" = "#E08214", "B" = "#E08214"))
Question: Hi there,

I have ran TopHat successfully with the fusions.out also. But I am unable to understand the directory structure i.e. how to create the directory structure although I have gone through the tutorial but always arising with errors in Tophat fusion.

The created directories are:

./blast_human/human_genomic*, nt* and other_genomic*
tophat_sample11 (containing the output files from TopHat2)
hg19
ensGene.txt
refGene.txt
mcl.txt
./blast/blastall and blastn

whenever I ran the command:

./tophat-fusion-post -o ./FUSIONOUT --num-fusion-reads 1 --num-fusion-pairs 2 --num-fusion-both 5 hg19

the error is:

[Mon Jan 21 16:00:55 2013] Beginning TopHat-Fusion post-processing run (v2.0.6)

[Mon Jan 21 16:00:55 2013] Extracting 23-mer around fusions and mapping them using Bowtie
samples updated

Traceback (most recent call last):
File "./tophat-fusion-post", line 2091, in ?
sys.exit(main())

File "./tophat-fusion-post", line 2059, in main
map_fusion_kmer(bwt_idx_prefix, params, sample_updated)

File "./tophat-fusion-post", line 315, in map_fusion_kmer
subprocess.call(cmd, stdout=open(output_dir + 'fusion_seq.bwtout', 'w'), stderr=open('/dev/null', 'w'))

File "/usr/lib64/python2.4/subprocess.py", line 412, in call
return Popen(*args, **kwargs).wait()

File "/usr/lib64/python2.4/subprocess.py", line 542, in __init__
errread, errwrite)

File "/usr/lib64/python2.4/subprocess.py", line 975, in _execute_child
raise child_exception

OSError: [Errno 2] No such file or directory

No log file is generated.

Thanking in advance!

Answer: The error was arising because Bowtie wasn't found in the path. So the path need to be defined as

PATH=$PATH:/path/to/bowtie

It worked for me...
Question: Dear Friends Hi (*'m not native in English!*)

I want to search for orthologs locally using EggNOG in my **fish** transcriptome samples.

First I have download [fiNOG][1] which is specially for fishes from [EggNOG][2],

then built a HMMER database using hmmpress (i.e. `cat fiNOG_hmm/*.hmm > fishDB.hmmer`), and run hmmpress fishDB.hmmer.

then I intend to run such script : `hmmscan --cpu 24 '/home/fiNOG_hmm/fishDB.hmmer' '/home/Transcriptome.fasta`'

**My question:** Can I use my transcriptome assembly .fasta file **directly** in this script or I must convert (translate) it into protein in advance (which tool is better for this job? *[Transdecoder][3] ?* )?

Thank you

[1]: http://eggnogdb.embl.de/download/eggnog_4.5/data/fiNOG/
[2]: http://eggnogdb.embl.de/#/app/seqscan
[3]: https://transdecoder.github.io/

Answer: There is now a better resource for functional annotation using eggNOG orthology. Check these links:

*Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper*

- http://eggnog-mapper.embl.de (online tool)

- http://biorxiv.org/content/early/2016/09/22/076331 (method description and benchmark)

- https://github.com/jhcepas/eggnog-mapper (eggnog-mapper tool. Use the option --translate if using nucleotide seqs)
Question: Hi everyone,

we have a new Biolinux system in the lab (Biolinux 8 which is Ubuntu Linux 14.04 LTS base) and I am struggling to install the program [Vicuna][1]. I am wondering if anyone could help me to sort it out. I think the major problem is the dependency from ncbi_tools++.

I have followed the installation for external users and downloaded [the NCBI toolkit][2], then configured it specifying the path, and finally `make` and `make install`.

The problem is that every time I am trying to run 'make' after the configuration of the ncbi_tools++ I get various make errors. And if I try to ignore these errors and install vicuna, the program will not run (I just gave it a try anyway).

I have tried to install the program both in bash and zsh environments as well as installing it in different locatoions without much luck.

I know Biolinux have already a version of the ncbi_tools

> ncbi-tools-x11 6.1.20120620 NCBI libraries for biology applications (X-based utilities)

But I am not sure if this is creating conflict with what I am trying to install.

If anyone could help me it would be much appreciated.

Any advice will be welcome

Regards

F.

[1]: http://www.broadinstitute.org/scientific-community/science/projects/viral-genomics/vicuna
[2]: ftp://ftp.ncbi.nih.gov/toolbox/ncbi_tools++/CURRENT

Answer: Let's put on the magic hat and see if we can guess what the problem is without seeing the errors...

If make works, but make install doesn't, then very likely it's a permissions problem.

The other common issue is missing headers.

That's about as far as the crystal ball with take us. If you post the errors from make, this should get a lot easier :)
Question: I want to count number of the variants called on each chromosome for each sample from a multi sample vcf file.
Any help would be really appropriated.
Thanks
Ram

Answer: Split your vcf file by sample and count how many times chromosom appear in each file .

FILE=yourfile.vcf
for sample in `bcftools query -l $FILE`
do
bcftools view -c1 -H -s $sample -o ${sample}.vcf $FILE
cat ${sample}.vcf |cut -f1|uniq -c > ${sample}.count
done
Question: I have few motifs whose location on a given DNA sequence has been identified using MATCH (from TRANSFAC). MATCH provides the location of these matrices in a numerical way as motif 1 236 (+) (for an example - motif 1 is present at location 236). I want to get a graphical output by feeding the length of DNA sequence and the the position of matrices so that its easy for presentation and comparing two different sequences with almost similar motifs

Thanks in advance for any suggestions

Answer: Have you tried the features and annotations options in Jalview
Question: Hi,

Is there a way to force GATK haplotypecaller to output all genomic positions in the vcf; not only variant ones ?

Example here :

Top : Raw bam file

Middle : -outbam using : `gatk HaplotypeCaller -R $genome -O test.vcf -I in.bam -ERC GVCF --max-alternate-alleles 3 -bamout test.bam -L interval.bed --output-mode EMIT_ALL_ACTIVE_SITES`

bottom : -outbam using `gatk HaplotypeCaller -R $genome -O test.vcf -I in.bam -ERC GVCF --max-alternate-alleles 3 -bamout test.bam -L interval.bed --output-mode EMIT_ALL_CONFIDENT_SITES`

You can see that the region without variant sites are not processed by haplotypecaller ; even with EMIT_ALL_ACTIVE_SITES and EMIT_ALL_CONFIDENT SITES.

Any idea how to tune haplotypecaller to output all sites with enough coverage ?

gatk version : 4.1.7.0

![enter image description here][1]

[1]: https://i.ibb.co/Dtk8WwK/igv.png

Thanks

Answer: `gatk HaplotypeCaller -R {reference} -I {bam_path} -L {bed_path} -ERC BP_RESOLUTION -O {out_vcf}`

Maybe look at BP_RESOLUTION
Question: FASTQC, for example, doesn't seem to have a publication associated with it. How would you cite it?

Answer: Just to say, thank you for wanting to cite tools. From our investigations, we've found that over 2/3 of the papers that mention Ensembl do not cite our papers. No idea how many more than that use our stuff and don't even bother to mention us. It's as if people think that building a bioinformatic database or tool is not real science, so doesn't need acknowledgement. In science, citations are currency, and using someone's work without citing them is essentially theft.
Question: I am trying to run the R command readVcf in R, it shows function not found.
I have already downloaded the package "VariantAnnotation". I don't know if it helps.
Does anyone have any idea?

Answer: Have you loaded it, `library('VariantAnnotation')`!!
Question: Hello everyone!

I would like to know if there is a way to get all the id's from BioSample. I already tried a link that i saw in another post that works for BioProject (https://www.ncbi.nlm.nih.gov/bioproject/browse/) but it doesn´t work for BioSample.

I also tried to download the summary of BioSample using this example from BioProject (ftp://ftp.ncbi.nlm.nih.gov/bioproject/summary.txt)

Also, can this be solved programmatically, i.e. using Eutils or EDirect?

What i really want to know is, simply, how many id's exist and how can i search for a list of them.

Thank you very much!

Answer: using my tool XsltStream http://lindenb.github.io/jvarkit/XsltStream.html and the NCBI Biosample XML dump

https://gist.github.com/lindenb/7d7397f6b11d140ed9c407e3e91a4c92
Question: Hi :)

A few weeks ago I was looking for a tool that would help me get "DNA composition" statistics for my sequencing. Something that would give me a dataset which which I could ask questions about GC bias, or over-represented sequences, motifs, etc. There are tools to answer each question specifically, but I was looking for something more general from which many analyses could be built on. This led me to k-mer counting, and all the down-stream tools which leverage k-mer result files.

k-mer counting tools are pretty cool, but all the ones I tried had some draw backs; like high RAM requirements, very long run-times (although the latests bloom-filter based tools seem to mitigate this somewhat), but most importantly requiring a specific k-mer size to be chosen. I really wanted all mers in the dataset so I could look at ''GC' and 'GCG' and 'GCCGACGGACGAC' without having to re-run any analyses. I couldn't find a tool like this after a brief search, so gave up and wrote my own in NumPy based off suffix-arrays.

Two weeks later I have a functional program in the sense that I get results, but before I invest any time making it usable for others, I thought I should investigate further if there are tools which already do this. Making a suffix array was a nice learning experience for me so I haven't lost anything if such tools already exist - and if they do i'd love to compare performance characteristics - but if not I might consider tidying up the code and making proper documentation. Does anyone know of such tools?

Thank you so much, and happy Diwali :)

Answer: For fixed k-mer, use kmc2 and DSK, especially when you have a fast disk. They are faster and use less memory than jellyfish.

If you don't want to use fixed k-mer, the right data structure is suffix array, FM-index or something equivalent. You can query any sequences up to the read length. Nonetheless, for a particular k-mer size, it is usually not as efficient as k-mer based methods.
Question: Dear lazyweb,

do you know any tool/script to visualize the different outputs of GATK DepthOfCoverage ?

P.

Answer: I wrote these scripts some time ago:

https://github.com/johandahlberg/Scripts/blob/master/plottingLociVsCoverage.R

https://github.com/johandahlberg/Scripts/blob/master/plottingCumulativeCoverage.R

R really isn't my strong suite, but at least its something, and they served my purpose at the time. I would be happy to accept any pull requests to improve on them, though.
Question: Hi
I am new in python and want to see how this can be done in python (I can do this in R). I have a text file `myfile.txt` with one column and thousands of rows as shown below. I want to convert this to fasta `result.fasta` format as shown below. How can I do this in python?

`myfile.txt`

ATGTGTGGTTTTCCCCC
ATTGGCGGGGTTTTTCAGGGG
ATGGGGGGGCCCCCCCCAAAAAA
TTGGTGGGGGGGGGGGGAA

`result.fasta`

>1
ATGTGTGGTTTTCCCCC
>2
ATTGGCGGGGTTTTTCAGGGG
>3
ATGGGGGGGCCCCCCCCAAAAAA
>4
TTGGTGGGGGGGGGGGGAA

Answer: #!/usr/bin/env python

n = 0
with open('myfile.txt', 'r') as f:
for line in f:
n += 1
print('>' + str(n) + '\n' + line.strip())
Question: Hi List,

Is there any way to get RNAfold predicted secondary structure in text format? It provides a PDF and PNG file but is it possible to get something like structures in this figure: [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3787264/figure/F3/][1]

Bade

[1]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3787264/figure/F3/

Answer: I had a Perl Module to output ASCII text of **miRNA precursor secondary structure**.

RNA::HairpinFigure - Draw hairpin-like text figure from RNA sequence and its secondary structure in dot-bracket notation.
http://search.cpan.org/dist/RNA-HairpinFigure/lib/RNA/HairpinFigure.pm#synopsis

>hsa-mir-92a-1 MI0000093 Homo sapiens miR-92a-1 stem-loop
CUUUCUACACAGGUUGGGAUCGGUUGCAAUGCUGUGUUUCUGUAUGGUAUUGCACUUGUCCCGGCCUGUUGAGUUUGG
..(((...((((((((((((.(((.(((((((((((......)))))))))))))).)))))))))))).))).....
---CU UAC C U UU
UUC ACAGGUUGGGAU GGU GCAAUGCUGUG U
||| |||||||||||| ||| |||||||||||
GAG UGUCCGGCCCUG UCA CGUUAUGGUAU G
GGUUU --U U - GU
Question: Hey guys,

I have an accession number, how I can convert the sra to fasta in Windows?

Answer: You can use [SRAtools][1] for this purpose. I think it is available for both Windows and Linux. Once installed, you can download/convert sra files to fastq as follows:

```
fastq-dump myfile.sra #for single end
fastq-dump --split-files myfile.sra #if paired end
```

[1]: http://www.ncbi.nlm.nih.gov/books/NBK158900/
Question: Go to ucsc genome browser, and enter the coordinate chr9:131,939,104-131,939,143 (hg19). It's part of the gene IER5, going from right to left. If you take a look at the DNA sequence CAG, you will see below that it codes for V (Valine). But, doesn't CAG code for Gln, and GAC code for Asp? Why does CAG correspond to Val on the browser?

Am I missing something terribly?

Answer: your gene is mapped on the reverse strand (3' -> 5') , the codon is not CAG but the complement: GTC (valine)
Question: Hi everyone; this is my first question on the forum.

How can I compare if two fasta files contain the same sequence headers?

Does any BioPython module exist for doing this?

Thanks in advance,
peixe

Answer: In BioPython you could do the following:

from Bio import SeqIO

def get_ids(fname):
reader = SeqIO.parse(fname, 'fasta')
ids = map(lambda x: x.id, reader)
return set(ids)

s1 = get_ids('f1.fasta')
s2 = get_ids('f2.fasta')

print 'Common:', s1 & s2
print 'In f1 and not f2:', s1 - s2
print 'In f2 and not f1:', s2 - s1

prints for my testcase:

Common: set(['a', 'c', 'b'])
In f1 and not f2: set(['x'])
In f2 and not f1: set(['y'])
Question: Hi!

Do you know how can I filter out supplementary alignments from a bam file? I was reviewing http://broadinstitute.github.io/picard/explain-flags.html and I am aware that the flag for this kind of alignments is "2048". However, depending on another features (ex: paired read, second in pair, etc), the flag can vary.

So, I am not sure about how can I filter out this kind of alignments.

Answer: `samtools view -F 2048 -bo filtered.bam original.bam`

You don't have to care about what other flags are set, the `-F` option will filter any entry with bit 2048.
Question: I want to compare data two different table like that using R.

Table 1
X1
a;b
c
d;e

Table 2
Y1
a uuu uu
b vvv vv
c xxx xx
d yyy yy
e zzz zz

Ouptput Table 3
X1 Y2
a;b uuu|uu<\br>vvv|vv
c xxx|xx
d;e yyy|yy<\br>zzz|zz

Answer: Here is the dplyr way. Make sure that you provide column names as given in `t1` and `t2`

library(tidyverse)
t1 <- tibble(X1 = c("a;b","c","d;e"))
t2 <- tibble(Y1 =letters[1:5] , Y2 = c("uuu","vvv","xxx","yyy","zzz") ,Y3 =c("uu","vv","xx","yy","zz"))

t3 <- t1 %>%
mutate(mm = map(X1, ~(strsplit(.,split = ";"))[[1]] %>% as_tibble)) %>%
mutate(nn = map(mm , function(.){
col_binds <- left_join(.,t2, by = c("value" = "Y1")) %>%
unite(col = "comb",sep = "|" , Y2:Y3)
row_binds <- col_binds %>%
summarise(value = paste0(value , collapse = ";") , times = paste0(comb , collapse = "<\br>"))
return(row_binds)
})) %>%
select(nn) %>%
unnest() %>%
dplyr::rename(X1 = value , Y2 = times)

## output
t3
# A tibble: 3 x 2
X1 Y2

1 a;b "uuu|uu<\br>vvv|vv"
2 c xxx|xx
3 d;e "yyy|yy<\br>zzz|zz"
Question: Hi!
I am trying to run the `drugInteractions` function in the MAFtools pipeline, but it seems that it is simply not there.
Here is the message I am getting:

> library("maftools")
> druggable = drugInteractions(maf = STES_maf, fontSize = 0.75)
Error in drugInteractions(maf = laml, fontSize = 0.75) :
could not find function "drugInteractions"

According to the authors, this function compiles info from the [DGI-db][1]. So I checked whether this funtion was borrowed from DGI, but I had no success.

Did anybody face the same issue? Would anybody know how to fix it? Not ideal, but recommending an alternative package for this kind of analysis would also help.

Thanks a lot in advance!

[1]: http://www.dgidb.org/

Answer: Problem solved!

There is no logic behind it, really. But it worked.
Well, simply get it directlly from GitHub:

devtools::install_github('PoisonAlien/maftools')

I have tried that on the day I posted the issue here, but didn't work. I did that again today, and it worked. :|

In any case, the author was very responsive! Very good feedback time. Thanks!
Question: Hi everybody!

Recently I became familiar with a new Gene Ontology plotting library in R called "GOplot". There is a tutorial in the developer site but from the very first point i have some problems with it. I put the tutorial in this question so you can review it: [http://wencke.github.io/][1]

1) There is no function in R like:

install.github('wencke/wencke.github.io')

2) If you skip this part you go through making dataframes. I don't know how to make mentioned dataframes in **toy example** part of tutorial for my microarray analysis. I would be really appreciate if you can help me with the second problem specially.

Thanks a lot.

[1]: http://wencke.github.io/

Answer: 1. You have a typo, try: `devtools::install_github("wencke/wencke.github.io")`, note that you need to have [devtools package][1] installed.

2. Data is within the package, you get them as noted in the tutorial `data(EC)`, it is a list of dataframes.

> **The toy example**
> GOplot comes with a manually compiled data set....

[1]: https://cran.r-project.org/package=devtools
Question: Hi everyone,
I have vcf file in dosage format like 0/0, 0/1 and etc. After filtering for quality, I need to impute the missing value using BEAGLE software so that my variants can be in float numeric format in order to combine with phenotypic data for further analysis. Apart from BEAGLE, if there is any other software that is much easier to use, I would not mind.
Could someone help me with script and guide me through imputation, please.

Answer: Which imputation program are you going to use?
Question: Dear All
I have a combined VCF file of few individuals. In this VCF instead of having normal CHR1, chr2 notions for chromosomes it is having the chromosome information as

gi|996703411|ref|NW_015379183.1|, gi|996703411|ref|NW_015379175.1

In this notion `NW_015379183.1` corresponds to a specific Chromosome. The same is true for its positions, If I have the chromosome numbers for all `gi|996703411|ref|NW_015379183.1|` sort of notions how I can replace the chromosome names to the original names.

Answer: use `bcftools annotate` https://samtools.github.io/bcftools/bcftools.html

with

> --rename-chrs file
> rename chromosomes according to the map in file, with "old_name new_name\n" pairs separated by whitespaces, each on a separate line.
Question: When I convert a SOAP Paired end output to a SAM file using the soap2sam.pl script(http://soap.genomics.org.cn/down/soap2sam.tar.gz) my SAM file headers turn out to be empyty.

Is there a way I can add headers to the SAM file ?

Answer: You can add the headers using samtools view:

samtools view -bT your_reference.fasta your_input.sam > your_input.bam

Use `samtools view -H your_input.bam` to check whether the header was added.
Question: My understanding of SAM files and the format is fairly good, but there are some things I haven't quite grasped. I'm not sure how obvious you all may find these questions, but they're what have come to mind.

I'm interested in recovering the original sequenced read after some alignment has been done. I'd like to know which pieces of the read/read segments I need. How do I know what the full read sequence was that came off the sequencer? Can I reconstruct it by connecting the read segments in the same template together? If so, what is the template? It's not the *read* is it? The SAM format doesn't suggest that it is; it says the template is some DNA/RNA fragment.

Here are some questions:

- "What is the difference between the read I get from sequencing and the read segments I see in a SAM file?"
- Intuition tells me that read segments are *mapped portions of the larger read*, but are they arbitrarily segmented in the SAM presentation?
- Are segments contiguous? Can they also be non-contiguous?
- Can I reconstruct the full read from the multiple read segments?
- How does *template* correspond to a sequence *read*?

I'm very grateful for any clarification I can get on these questions.

Answer: The given read from your query file (fastq file) can match to multiple locations in genome. You can check this with NH flag in sam file. The reads can also overlap each other as they are sequenced from DNA fragments and you can find this by comparing the mapping co-ordinates in sam file.

The aligner take only read sequence from fastq file for mapping to reference sequence. If you want use contigous sequence (contig), you need to assemble it first and then map with reference sequence. As the contig will be longer in length, you need to be cautious while using aligner.