cannin/biostars_qa huggingface dataset:
Question: Hi,

I get an error while making txdb from biomart for S. pombe. Can someone tell what am I doing wrong? It works when I give biomart="fungal_mart" . I am not sure if the txdb will be different?! Thanks . . .

txdb<-makeTranscriptDbFromBiomart(biomart ="fungi_mart_23" ,dataset="spombe_eg_gene" ,host="fungi.ensembl.org")
Error in useMart(biomart = biomart, dataset = dataset, host = host, port = port) :
Incorrect BioMart name, use the listMarts function to see which BioMart databases are available

Answer: The error message tells you the problem. The mart *fungi_mart_23* is not associated with the host *fungi.ensembl.org*.

```
listMarts(host = "fungi.ensembl.org")

biomart version
1 fungal_mart Fungal Mart
2 fungal_variations Fungal Variation Mart
```

If you do not specify a host, the default is *www.biomart.org*, which does contain *fungi_mart_23*.

```
mart <- listMarts()
mart[grep("fungi", mart$biomart),]

biomart version
5 fungi_mart_23 ENSEMBL FUNGI 23 (EBI UK)
6 fungi_variations_23 ENSEMBL FUNGI VARIATION 23 (EBI UK)
```
Question: Hello,

I used Cufflinks (v.2.0.2) to assemble a complete transcriptome without providing any GTF references just to check the quality of assembly and it seemed to work so well, but I am worried several terminal exons were not entirely reconstructed (see attached figures in which I compared my cufflinks results with RefSeq annotations using IGV browser) although read coverage was quite good. I tried to change several parameters (disabling low-coverage 3' trimming; decreasing minimum reads per transfrags, min-isoform-fraction, or small-anchor-fraction) but the problem persisted.

Did somebody have the same problem or know how to fix it?

Thanks,
Jorge

Answer: Hi everyone,

I realized that it is a problem of Cufflinks version (v.2.0.2). Using newer versions, terminal exons are correctly assembled.
Question: as title, there is no overlap region in bed file. Thanks

Answer: You can do it with the following command line:

cat file.bed | awk -F'\t' 'BEGIN{SUM=0}{ SUM+=$3-$2 }END{print SUM}'
Question: Hi !
I'm currently a student and I have a hard time understanding some basics of bioinformatics; I'm currently learning about *alignment*, *filtering*, *variant calling* and such so my question might look silly but here it is anyway.

I have some trouble about how you work with paired-end sequencing files and what does it means to be paired-end.

After taking a look on the Internet I found an explanation of what is paired-end sequencing (tell me if I got it right):

*For what I understood, a paired-end sequencing is just done by sequencing from A to Z and then from Z to A. Which will provide two distinct datasets, one for each direction.*

My question is, when you are doing some alignment with tools like *BWA*, *TopHat* or whatever, do you have to reverse one of the two dataset or not ? Because, for instance, If I wanted to find a consensus sequence (or the position specific score matrices), if half of the data are in the wrong direction wouldn't it be completely wrong ?

Completely unrelated: I've also heard that TopHat should be used over BWA for aligning RNA, do you know why ?

Answer: Always easier to illustrate with an image, from [here][1]. The grey represents the fragment, and each end of the fragment is sequenced. This allows more accurate mapping, particularly of repetitive regions. There's also a great animation [here][2] that illustrates the concept of Illumina paired end sequencing. As @h.mon stated, most programs will have parameters to deal with paired-end sequencing, and seriously, stay away from Tophat. [STAR][3] or [HISAT2][4] are current alternatives

![enter image description here][5]

[1]: https://www.illumina.com/science/technology/next-generation-sequencing/paired-end-vs-single-read-sequencing.html
[2]: https://www.youtube.com/watch?annotation_id=annotation_228575861&feature=iv&src_vid=womKfikWlxM&v=fCd6B5HRaZ8
[3]: https://github.com/alexdobin/STAR
[4]: https://ccb.jhu.edu/software/hisat2/index.shtml
[5]: https://assets.illumina.com/content/dam/illumina-marketing/images/science/v2/web-graphic/paired-end-vs-single-read-seq-web-graphic.jpg
Question: Hi,

I have a bam files produced from mapping by LifeScope Tool (mapping of SOLiD reads). I would like to know the length of the reads that were mapped to the reference genome based on the information present in bam files. Is there a tool which can give me the stats of such things from bam files?

Thanks!

Answer: This will take the first million reads and give you counts of each read length. To look at more of the bam, get rid of the "head" command:

samtools view file.bam | head -n 1000000 | cut -f 10 | perl -ne 'chomp;print length($_) . "\n"' | sort | uniq -c
Question: Hi,

I have around 20 BAM files and each of them have proper bam header including RG ID, Sample, Library, Platform unit, Description and Platform. I was reading the GATK manual and just realized that the only accepted names for platforms are 454, LS454, Illumina, Solid, ABI_Solid, and CG. But I wanted to be quiet specific when I was converting SAM to BAM and used "Solid5500XL" as the platform name.

I know there are AddorReplaceReadgroups and ReplaceSamHeader coomand in Picard which I can use to edit the BAM header. But it will take some time as I have many large files. Is there an easy way to do so without rewriting the BAM files.

Thanks in advance

Answer: I think you can do this with [samtools][1], first extract the header into a text file, then correct the value and then change the BAM header. Something like this:

```
for BAM in *.bam
do
samtools view -H $BAM > header.sam
sed "s/Solid5500XL/Solid/" header.sam > header_corrected.sam
samtools reheader header_corrected.sam $BAM
done
```

I haven't tested it but I think it must work for you. Please let me know if that works (who knows, maybe some day I'll need it).

[1]: http://samtools.sourceforge.net/
Question: I have ~70k dbSNP ids I need to get Exac minor allele counts and minor allele frequencies for.

Supposedly Exac provides this information for variants, but I can't find it, either in the JSON results obtainable through http://exac.hms.harvard.edu/rest or through Exac's web interface.

(BTW, the information needs to be from Exac, to be procedurally consistent with other annotations I'm working with.)

Answer: Hello,

you can use the [dbsnp endpoint][1] from Exac's REST-API, or the
[variation endpoint][2] from ensembl's REST-API. In both cases you have
to parse the json for the allele frequencys to find out which allele is
the minor one and what it's value.

fin swimmer

[1]: http://exac.hms.harvard.edu/#rest-dbsnp
[2]: http://rest.ensembl.org/documentation/info/variation_id
Question: I need to find the nucleotide diversity at each position of an alignment. Is there a way to do this with available software ? I would also be comfortable if there is a way to do this using R

Answer: You can use WebLogo. The output is either an image (in several formats) or a text file with relative frequency of each nucleotide/amino acid at every position in the alignment.
Question: When learning bioinformatics, I'm constantly coming across new terms which I need to find the meaning of. Currently, to learn the meaning of terms, I Google the term, and usually read a Wikipedia article.

The question: Is there a curated, high-quality dictionary of biological terms available that beats random Google searches in terms of quality/information content?

Answer: How about a collection of ontologies? Check out obofoundry. That is maybe as close as it gets at the moment, use the search to find terms. Also machine readable, useful for data integration and automatic annotation.
Question: I am looking for a single file that lists species name for every taxonomy id in NCBI. Is there such a file? I was unable to find a flat file in the taxonomy ftp site. I have been using the site below to manually do the matchings so far but it is quite un-scalable

http://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi

Answer: I think you can find that information in a file called "names.dmp" that you can derive from ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz. Take a look at this readme: ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump_readme.txt
Question: Aloha,

I am having trouble figuring out how to remove everything after the last '_' in the sequence headers of a fasta file.

I would like this following series of headers

>ART01B_100_M7_ID100005_1
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325_189
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005_46
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT

to look like this:

>ART01B_100_M7_ID100005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT

because there are various '_' in the sequence headers, this command, for example, won't work:

cat Alt_MACSE_Output.fasta | awk -F _ '/^>/ { print $1"_"$2"_"$3"_"$4 } /^[A-Z]/ {print $1}' > Alt.fasta

Can anyone help me please?

Answer: awk 'BEGIN{OFS=FS="_"}{if(/^>/){NF--}}{print $0}' in.fa > out.fa
Question: Hi,
please, how can I visualize nucelotides for every position. I have a percentage of bases but now I don't know how to create a plot.

fastq <- readDNAStringSet("https://d28rh4a8wq0iu5.cloudfront.net/ads1/data/ERR037900_1.first1000.fastq","fastq")
fastq
freq <- alphabetFrequency(fastq, as.prob = T,baseOnly=T)

I would like something like this (sorry for quality):

![][1]

[1]: https://image.ibb.co/cQPULG/base.png

Answer: You need to change your code. It is not alphabet frequency. It is consensus matrix. Example code and image below:

library(Biostrings)
fastq <- readDNAStringSet("sample.fastq","fastq")
## At each position, base frequency
afmc=consensusMatrix(fastq, baseOnly=T,as.prob = T)
tafmc=t(afmc)
matplot(tafmc[,-5], type="l", lwd=2, xlab="Read Length", ylab= "Base frequency at each position")
legend(legend = colnames(tafmc)[-5],"topright",col=1:4, lty=1:4, lwd=2)

5th column is others (other than A,T, G and C).
Question: Hi, I have a suspicion one or more of my fastq files is corrupted: some reads have sequence and quality of different lengths, for example. Would this python script work to detect a problem?

import sys
from Bio import SeqIO
for record in SeqIO.parse(sys.stdin, "fastq"):
pass

And I use it like this:

pigz -dc -p 3 serum-3_2.fq.gz | python check.py

Will it raise an error if one of the reads is non-valid? Or is there a better quick way to check this?

Answer: Yes, it raises an error if sequence and quality strings have different length.
Question: I wanted to see a particular sequence in our re-sequencing data using IGV. However the .bam file is too big (10Gb) to handle. Since I am interested in only a small region in the genome, I wondering if there is any way to generate a .sam/.bam file only for that region. Thank you very much in advance!

Chenglin

Answer: This has been answered many times before. You can specify a small region of your interest and create a small bam file as shown.

samtools view -bh chr1:100-200 > small.bam
Question: Hi, I'm a computer scientist working on PPI networks. I'm trying to understand the concept of synonym. How are synonyms defined and how can one conclude that 2 gene names are synonyms?

I thought they would be the same or almost same, having very similar sequence and interactions but I found counter examples. For example I'm wondering how this situation does exist:

I found a gene in Homo Sapiens called MEN1 (official symbol). When I look at its synonyms, I see SCG2 which is the official symbol of a gene also in Homo Sapiens.

There are 2 genes. They are both in Homo Sapiens and they are synonyms. Their interactions are completely different.

I could not find much information about this issue. I'd like to know the definition of synonym and if I can assume two genes that are synonyms are the same for some specific cases.

Answer: Synonyms are a curation issue. Alternative names are generated from a variety of situations: some are different labs stumbling on to the same gene, some are simple differences between characters (MEN1, MEN 1, MEN-1). Sometimes large-scale projects generated some form of coded designation (such as 1200015E08Rik, a mouse example).

Alternate names are collected by curation teams at various places. Some come from papers read by curators. Some would have come from computational sources. Most of them are shared around, but there may be some sources that have some name that others don't.

There are nomenclature standards that are supposed to be used, and there are committees that establish the rules and assign official names and symbols. However, some people refuse to use those names in publications and revert to their preferred name.

The official name/symbols shouldn't have duplicates. But the synonyms might--so what you are seeing with SCG2 could certainly be 2 different genes.

EDIT: If you are new to this arena maybe you should see some nomenclature stuff. Here's the mouse stuff, maintained by Jax, where I learned about it:
http://www.informatics.jax.org/mgihome/nomen/index.shtml

Here's the HUGO site for human http://www.genenames.org/
Question: I was wondering whether anyone had trouble with using RECON when running REPEATMODELER.

I have downloaded the latest version of RECON (1.08). When prompted by ./configure in REPEATMODELER, I give the file path of RECON as /$RECON/bin. However, the message "RECON-1.08/bin does not appear to contain a complete set of files". The /bin file is empty, however, do I have to put executable files in it - if so, how do I do this?

Does anyone know how I solve this issue?

Thank you

Answer: Have you compiled RECON?

If not, open a terminal in the src folder, then run " make" and then "make install".
Next in the script "recon.pl" in the scripts directory -- on the third line, add
the path to the binaries (the bin directory here) between the double
quotes.

I know this answer is too late for you, but it might help others.
Cheers.
Question: I have a biomart xml query and I am going to obtain the data with curl and the RESTfull de biomart.

In the dataset it is defined the specie and ensembl but not the version of ensembl to use.​

Question: where I specify the build?

Answer: I think you have to select the correct host for your data: http://grch37.ensembl.org/biomart/martview
Question: Can someone compare and contrast the following scientific workflow frameworks?

Taverna http://www.taverna.org.uk/

Pegasus http://pegasus.isi.edu/

Answer: What sort of criteria did you have in mind? Yes, they're both (fine) workflow systems but they're rather different from each other.

Pegasus is focused on coordinating programs that transform files into other files, which is definitely a powerful paradigm. In tooling terms, it works through command line tools only: if you're working with complex workflows, you'll _definitely_ want to have an extra tool around it to generate the workflow for you from some higher-level description. The tooling is such that it is relatively easy to use it behind websites, with other programs, etc. Because you typically generate the workflow programatically, you can easily adapt the workflow according to your particular requirements (provided you're a programmer).

Taverna is much more of an integrated system, and it works at a higher level of abstraction — not files, but rather the data that you might find in them, and not programs, but rather the operations you might do with them. It comes as a GUI workbench for creating the workflows, running them and visualizing them, a command line tool for just running them, and a server for running workflows when you need more power (e.g., to support a portal). The workflows are comparatively awkward to edit by mechanisms other than using the workbench.

To look at the other things you mention in your comment:

workflow and makefile concept features — Neither system is make and neither works like it. Yes, they all support branching, parallel execution and data merging, but some of the other things you mention under this just don't make sense.

ease of development — It really depends on who you are and and how you like to develop. In my opinion, Pegasus is perhaps better for workflow creators who are principally programmers (because it has a less complex workflow format to generate) and Taverna is perhaps better for workflow creators who are principally scientists (because it comes with a GUI tool for the job).

support — Both are supported. Both are producing new versions. There's multiple ways to get support in both cases. I'd call this one a score draw…

cluster support and tracking — I'm not quite sure what you mean by this. The problem is that when you're dealing with workflows that may potentially be accessing many remote services (that might themselves be cloud-hosted), the meaning of clustering gets very… slippery… and I really don't know what you mean by tracking.

workflow visualization — Taverna does this; has done for many years. There may be a separate tool for this with Pegasus, but I don't know it (and the different pattern of workflow that is often associated with Pegasus is considerably more difficult to visualize for various reasons, notably that Pegasus workflows are frequently using very wide span-outs; this is a problem that Taverna avoids by representing such things as a single processing element that has multiple data items flowing through it).

_Declaration: I'm one of the Taverna developers._
Question: First consider the naive non-automatic approach to this question.

Lets use Escherichia coli str. K-12 substr. MG1655 as an example it has the taxonomic ID 511145. In the gnome DB of NCBI I can find this entry:

https://www.ncbi.nlm.nih.gov/genome/?term=txid511145[Organism:exp]

The first line says -> Reference genome: Escherichia coli str. K-12 substr. MG1655. Which links to:

https://www.ncbi.nlm.nih.gov/genome/167?genome_assembly_id=161521

Which is what I want. I can no use the RefSeq Accession to get the genome.

Now lets try this automatically:

esearch -db genome -query "txid511145[orgn]" | esummary

The result is the following:

167
Escherichia coli
Bacteria

Gammaproteobacteria
A well-studied enteric bacterium
12319
PRJNA12319
Complete
3
10
0
ASM2632v1
GCA_000026325.2
1181661
1998/10/13 00:00

900
484
9772
0
562

The first problem is that this querry brings me to species level (see TaxId 562) and not to the far more specific level of my query (strain).

The specific question in this case would be how can you link the TaxId 511145 to the genome_assembly_id 161521 using an esearch querry. Obvious the connection is there as I can make the link on the website.

Answer: #!/bin/bash
TAXID="$1"
if [[ -z "$TAXID" ]]; then
echo "Usage: retrieveReference.sh NcbiTaxId"
exit
else
if [[ ! -e assembly_summary_genbank.txt ]]; then
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/assembly_summary_genbank.txt
fi
awk -v TAXID="$TAXID" \
'BEGIN{FS="\t"}{if($5=="reference genome" && $6==TAXID){print $20}}' \
assembly_summary_genbank.txt \
| awk 'BEGIN{OFS=FS="/"}{print $0,$NF"_genomic.fna.gz"}' \
| xargs -n1 wget
fi
Question: Hello,

I have what might be a simple question but I would really like to get this right and despite all my readings it is still not clear to me what the best option is.

I am using DESeq2 to analyse an RNAseq experiment. The design of the experiment is as follow (each sample is in triplicate):

- cell line A - untrt
- cell line A - trt
- cell line B - untrt
- cell line B - trt

The goal is to peform the following contrats:

- cell line A - untrt **versus** cell line B - untrt
- cell line A - trt **versus** cell line B - trt
- cell line A - untrt **versus** cell line A - trt
- cell line B - untrt **versus** cell line B - trt

My question is, in order to do things "correctly", should I create a separate DESeq object for each contrast using only the samples I will be using for each contrast, or should I have just a single object from which I do the contrasts that I need?

For your information, the difference between these cell lines is that "cell line A" is overexpressing gene X, while "cell line B" is overexpressing genes X and Y.

Any help would be very much appreciated...

Answer: "should I create a separate DESeq object for each contrast using only the samples"

Create one object.

and then for comparison,

do this -
res, comparison1<-results(dds, contrast = c("A.untreated", "B.treated"))
and so on.

You can list the contrast names in the object by doing resultNames(object).
Question: I'm new to Galaxy and bioinformatics tool.

I have Fastq data of chicken heavy and light chains and I want to do pca plot on them. I know I have to convert them to BAM file first but I don't know which tool to use - FastqtoSam or Bowtie2. One gave me unaligned BAM and another gave aligned BAM. Then, I used multiBAMsummary tool to see which file was correct but both BAM file gave me an error saying file '0.bam' does not have BAM or CRAM format.

Where did I go wrong?

Answer: That means it did not run. You can align the fastq files using `bowtie2`:

bowtie2 -x index -1 R1.fastq -2 R2.fastq -S bowtie.sam

Then convert the `sam` to `bam` using `samtools`:

samtools view -bS bowtie.sam -o bowtie.bam

Then `sort` it:

samtools sort bowtie.bam -o bowtie.sorted.bam

You can `index` the `.sorted.bam` by using:

samtools index bowtie.sorted.bam
Question: I'm trying to figure out what tools most people use for differential chromatin detection among samples (condition1 vs condition2).
I normally use diffReps for differential detection for ChIP-seq but haven't done it for ATAC. Any suggestions?

Answer: 1. Call the peaks, by merging all the replicates (bam files), with very less stringent criteria ( p-value of 0.1 ) such that you end up with union of regions ( peaks ) across two biological conditions, that have some signal enrichment.

2. Now either use the entire peak or center the peaks on summit with +/- 200bp and count the number of reads with in each peak from all replicates independently using featureCounts. ( You can easily covert the above `bed` file to `saf` format )

awk -v OFS = "\t" 'BEGIN { print "PeakID", "Chr", "Start", "End", "Strand" } { print "Peak_"NR, $1, $2, $3, "." }' test.bed > test.saf

3. Feed this matrix to DESeq2 or edgeR and get differential peaks. **Do some filtering of peaks like removing peaks that have very few counts across replicates, otherwise this will be a problem for multiple testing correction to include so many peaks .**

**Edit:**

Most of the differential binding tools run either DESeq or edgeR in the back end and they are suitable for TFs, as the signal is sharp with out any noise.

But with ATAC, the peaks have both characters of TF and histone modifications.

I would

1. Run macs2 with p-value of 0.1

2. Use the narrowPeaks file, either the entire peak (Better)or take +/- 200bp from the summit (only if peaks called with --summit option in MACS, so that broad peaks are split into multiple sub-peaks).

3. Quantify peaks using featureCounts.

4. Feed the matrix to DESeq2.

**You should also upload this file ( from step 2 ) to browser to get a feel of how the regions that you are using for DE analysis looks like. Most likely they will be okay.**