cannin/biostars_qa huggingface dataset:
Question: Hi,

I am trying to blast contigs to the C. elegans reference genome sequence with the following 2 commands:

makeblastdb -in celegans.fa -input_type fasta -dbtype nucl -title celRef -out celegans.db -taxid 6239

blastn -subject celegans.db -query celegans.contigs.fasta -out bristol.txt -outfmt 6 -html -num_threads 8

More details here:

makeblastdb -in output/genome/celegans/ref/seq/celegans.fa -input_type fasta -dbtype nucl -title celRef -out output/genome/celegans/ref/database/celegans.db -taxid 6239

Building a new DB, current time: 08/30/2017 16:35:01
New DB name: /Users/cr517/Documents/phd/project/output/genome/celegans/ref/database/celegans.db
New DB title: celRef
Sequence type: Nucleotide
Deleted existing Nucleotide BLAST database named /Users/cr517/Documents/phd/project/output/genome/celegans/ref/database/celegans.db
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 7 sequences in 1.03467 seconds.

blastn -subject output/genome/celegans/ref/database/celegans.db -query output/assembly/celegans/canu/bristol/celegans.contigs.fasta -out output/alignment/celegans/ref/pacbio/contigs/canu/blastn/bristol.txt -outfmt 6 -html -num_threads 8
Command line argument error: Argument "subject". File is not accessible: `output/genome/celegans/ref/database/celegans.db'

But when I check with ls, I do have database files in my directory:

ls output/genome/celegans/ref/database
celegans.db.nhr celegans.db.nin celegans.db.nsq

I do not understand why the database is not accessible. Your help would be appreciated.

Thanks.

Best, C.

Answer: You want to use `-db /path/to/celegans.db`. `-subject` is used for blasting two sequences against each other - I think both in fasta format.
Question: Hi,

I actually try to run ChromHMM on a new cell line for which I have the following markers (H3K4me1, H3K4me3, H3K9me3, H3K27ac, H3K27me3, H3K36me3, Dnase Hotspot and input) and I have some difficulties to retrieve right annotation corresponding to different chromatin state. Does somebody have some advices/tips to share?

Does somebody already try this software on a new cell line?

Thanks in advance,

Answer: It's an art more than a science. It takes a bit of staring at the emission probabilities. Use the OverlapEnrichment tool with a bunch of BED files of factors whose biological function is well-studied to help interpret the models. For example in your case high H3K4me1, low H3K4me3, and high H3K27ac for some state N suggests a label of "enhancer-like" for state N. You could use OverlapEnrichment with, say, P300 peaks (something not used in the model) to help justify that label. Enhancers, active promoters, and repressive chromatin are usually pretty straightforward to label. The other states that don't have an obvious succinct label get tricky and end up with awkward names. That's biology for you though -- defying simplicity!
Question: I have a VCF file, lets call it file A, which was created by performing variant calling using GATK on whole genome sequencing data. For SNPs that do no appear in the VCF, can I assume that the SNP is monomorphic for the reference allele (i.e all 0 encoding). Note that sites that are monomorphic for the alternate allele appear in the file. I ask because I need to merge this file (A) with another VCF (file B) and I'm not sure how to handle the variants that appear in B but not A. Can I just fill in the genotypes with all 0 for these variants that do no appear in A, or do I have to impute first to account for the possibility of low coverage in the region? Would be great if answers provide a external source as well so I have to point of reference when I consult my advisor because he thinks that sequencing data should never require imputation. Thanks!

Answer: I wrote a tool to fix this problem : http://lindenb.github.io/jvarkit/FixVcfMissingGenotypes.html

(but it is slow)

> After a VCF-merge, read a VCF, look back at some BAMS to tells if the missing genotypes were homozygotes-ref or not-called. If the number of reads is greater than min.depth, then the missing genotype is said hom-ref.

see also: https://www.biostars.org/p/119007 https://www.biostars.org/p/263309 https://www.biostars.org/p/276811
Question: In R, how to get the gene symbols of a GPL (like [GPL570][1])? There is no gene symbols in the annotation for some GPL from GEO. I have lots of GPL from multi platforms. for these with no gene symbol, they have `GB_ACC`, so more specifically, how to map `GB_ACC` to gene symbol in R?

Tools like [ailun][2] can do this online or WSDL or Perl. But not in R or other language I master. Are there any other tools?

[1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL5717
[2]: http://ailun.stanford.edu/platformAnnotation.php

Answer: Each GPL is provided by the submitter, so there is no "standard" approach for looking up gene symbols. You'll need to examine all the GPLs and determine what column(s) could be used for lookup.

Just a note that Bioconductor can be really useful for these mapping tasks and several (mainly Affymetrix) array platforms have corresponding Bioconductor data packages that can be used to facilitate ID mapping. A mapping from GPL to Bioconductor annotation packages is available via the GEOmetadb Bioconductor package or at this gist:

https://gist.github.com/seandavi/bc6b1b82dc65c47510c7
Question: I have been using `heatmap.2` for a while, but just discovered `pheatmap`. In `heatmap.2`, you can specify clustering settings via `distfun` and `hclustfun`. In `pheatmap`, you have `clustering_distance_rows` and `clustering_method`. However, if I set those parameters to use the same algorithms, the resulting heatmaps do not look similar. How can that be? Does `pheatmap` perform additional manipulations that `heatmap.2` does not?

My code:

```r
# pheatmap
pheatmap(vals, color=colors, scale="row", cluster_rows=T, cluster_cols=T, clustering_distance_rows = "euclidean", clustering_distance_cols = "euclidean", clustering_method = "complete")

# heatmap.2
hclust_fun = function(x) hclust(x, method="complete")
dist_fun = function(x) dist(x, method="euclidean")
heatmap.2( as.matrix(vals), scale="row", trace="none", dendrogram="both", Rowv=TRUE, Colv=TRUE, distfun=dist_fun, hclustfun=hclust_fun, col=colors)
```

Answer: Basically when you show scaled data, heatmap.2 scale data after clustering , whereas pheatmap scales data before clustering. I am guessing that makes the difference in the final output sometimes.
Question: I am building a large database for my analysis tool ~3 Terabytes.

The issue is distributing this database for people to use it.

A solution we have learned from working with sequencing cores is to ship a hard drive by FedEx. A possibly expensive and annoying proposition.

Another solution is to host the files for download, but this seems possibly expensive or possibly impractical for the audience. (On this point of practicality -- most Universities have gigabit level bandwidth, perhaps 3 Terabytes is no longer a ludicrous download size, simply a big one.) At 1 gigabit, I estimate the download would take 400 minutes.

Since this is a niche' problem, I am hoping someone may have had a similar issue and can provide advice?

Thank you
Jeremy Cox

Answer: You have already covered all bases. There is no reduced cost alternative here.

Amazon or Google storage (perhaps the cheapest data buckets if no compute needs to be done on the data). Storing it in one place and sharing as needed is the best solution. Google storage has [good access control options][1] that work in google compute.

[1]: https://cloud.google.com/storage/docs/access-control/
Question: Hi,

I'm currently using BamTools to parse a (paired-end) RNA-seq bam file for detecting possible chimera.

I'm focusing on reads where both pairs belong to different genes but I have a hard time to understand some of the BamTools::BamAlignment members.

For example, it is easy to get the starting position of the mate read using the "al.MateRefId" for the chromosome and "al.MatePosition" for the coordinate but I still lack its end position. For standard reads (both pairs belonging to the same gene) I can get the mate end-position using "al.Position + al.InsertSize" but is it still true for chimeric reads? If not, how is InsertSize computed and how can I get this mate end-position?

Many Thanks in advance

Answer: The insert size it specified by the aligner and for RNAseq where the paired-ends don't overlap or when they span a splicing event or other type of junction (e.g., fusion bounds), these may be unreliable. The end position of the mate isn't stored in the BAM format, so you need to find the mate and then get its end. I'm not familiar with the Bamtools API, so I can't say whether it provides a convenient method for this or not, though I wouldn't be surprised if it doesn't.
Question: I have a very big compressed vcf file with over 8000 samples. I wanted to extract out data for those samples whose name start with first three characters "TMS" so that I can have a new vcf file with those samples and their variants data only. Could anyone help me out with the vcftools or bcftools command to accomplish this?
Thanks

Answer: gunzip -c in.vcf.gz | grep "#CHROM" -m 1 | cut -f 10- | tr "\t" "\n" | grep ^TMS > samples.txt

bcftools view --samples-file samples.txt in.vcf.gz
Question: Hi,

I aligned a few samples using STAR to the genome provided in the Illumina iGenomes UCSC hg19 bundle ([here][1]) -- I used the provided gene feature (gtf2) file as is. Now, my motive is to calculate the gene and isoform expression levels using bedtools multicov (at the same time).

Use of the gtf2 file produces a file containing read counts per exon. I wish to compute gene and isoform read counts too, so I converted the gtf2 file to a gff3 file using using gtf2gff3 script from SO/GAL ([here][2]).

My first question is: Is it OK if the alignment is performed with gtf2 file but counted for reads using the gff3 file, keeping in mind that the gff3 file was converted from the gtf2 file?

My second question follows

I have read both these resources ([here][3] and [here][4]) but do not understand the differences between:

- exon vs CDS
- transcript vs mRNA

I know that with the process I described, it is possible to retrieve gene read count by selecting only the lines where feature=gene from the bedtools multicov output. What must I do for isoforms? I am confused by the semantics.

Thanks ahead of time and let me know if my post was not clear enough.

[1]: http://cufflinks.cbcb.umd.edu/igenomes.html
[2]: http://www.sequenceontology.org/software/GAL.html
[3]: http://mblab.wustl.edu/GTF2.html
[4]: http://www.sequenceontology.org/gff3.shtml

Answer: You were almost at the right place, the primary purpose of the website that hosts the GFF3 definition is to define the terminology (ontology) for each term. If you go to the browser and enter *exon* and *CDS* you will get the exact specification for each.

http://www.sequenceontology.org/browser/obob.cgi

Exon: *A region of the transcript sequence within a gene which is not removed from the primary RNA transcript by RNA splicing.*

CDS: *A contiguous sequence which begins with, and includes, a start codon and ends with, and includes, a stop codon.*

As for conversions, the GTF2 and GTF3 files differ only in the information contained in the 9th column.
Question: I have a read where sam flag is 69 (meaning Pair End Read, Unmapped and First in Pair) , Mapping quality is 0 and Cigar string is `*`. The sequence is not showed in the viewer (samtool tiview,IGV).

I'm wondering why a read is said unmapped whereas a position is given in the sam file. So the read can be mapped? When I tried to align with blat, I found one unique position (the same as in sam file).

So I don't understand ...Why this read is said unmapped whereas a position is given(chr12 112036669)? How to retrieve/use this read in my bam to display in viewer as a "good sequence"? Why is the read filtered? Because of the mapping quality? and why CIGAR string is set to `*`? Can I use it of one way or an another or is it a better idea to let this read stay in "garbage"? Thanks.

M02792:15:000000000-A98VA:1:1103:24705:18436 69 chr12 112036669 0 * = 112036669 0 GTCGCTGAAGCCCCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCCGCCGCCCGCGG CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG8FGGG@FF:F@FE=: RG:Z:A421

Answer: If an alignment is marked as unmapped then any position information is unreliable, since it may simply be a hold-over from an internal method that didn't get cleared. In the case of paired-end reads like this, typically an unmapped mate is given the coordinates of its mapped mate.
Question: Hello,

I am new to these technologies, and would like some more information to distinguish each. In what ways can CAGE and RNA-seq data complement each other? What is unique for each of the technologies?

For instance, given I have some data that supposedly also has matching CAGE data that is mapped and clustered. If these samples match the RNA-Seq samples, how do I know if the CAGE samples can complement the RNA-Seq analysis?

Thanks for any help!

Answer: In CAGE, only mRNA with 5' caps are captured and sequenced.

Usually there will be multiple TSS for a gene and some genes may have multiple promoters. In RNA-Seq you will get to know which genes/transcripts are expressed, because the mRNA is randomly fragmented and sequenced, but that does not tell the which promoter being used. CAGE-Seq on the other hand, sequences the 5' regions of mRNA, which tells you the precise location of TSS at 1bp resolution, which explains you about the promoter usage and promoter architecture.

In short, if you want to study the promoter architecture, you need to use CAGE "like" data. If you are only interested in expression, RNA-seq would be enough. If you have deep enough the CAGE data, you could study both the promoter architecture and the gene expression.
Question: Going through my RNA-seq samples, I've noticed this weird pattern for 'per base sequence content' graph.

[http://imgur.com/bhYIVJT][1]

There is a lot of noise in the first few base pairs and I was curious what it is and if I should be concerned for not reducing the noise before alignment.

[1]: http://imgur.com/bhYIVJT

Answer: Here is a [link][1] (for a post by Dr. Simon Andrews, author of FastQC) to put this topic to rest. Short answer: Nothing to worry about or mitigate.

[1]: https://sequencing.qcfail.com/articles/positional-sequence-bias-in-random-primed-libraries/
Question: Is there any tools or commands to sort bed format file? The correct order is chr1, chr2, chr3, chr4,..., chr22, chrX, chrY.

Answer: GNU sort can now sort in alpha-numeric order. See option `-V, --version-sort` of the following version "sort (GNU coreutils) 8.17".

Here is an example:

$ echo -e "chr10\t1\t2\tA\nchr2\t1\t2\tB\nchr1\t1\t2\tC"
chr10 1 2 A
chr2 1 2 B
chr1 1 2 C

And the result is:

$ echo -e "chr10\t1\t2\tA\nchr2\t1\t2\tB\nchr1\t1\t2\tC" | sort -k1,1V
chr1 1 2 C
chr2 1 2 B
chr10 1 2 A
Question: I have a nextflow workflow for which outputs are saved to a `results` directory. The nextflow `main.nf` script and `rscript.Rmd` are in the `workflow` directory.

When I run rscript render, the `rscript.Rmd` automatically generates the output in the `workflow` directory rather than the `results` directory as desired.

To remedy this, I added the parameter output_dir to the render function. However the output ended up at this path: `workflow/work/82/results/run_quality_report.html`. And there was error `Missing output file(s)`

publishDir params.outdir, mode:'copy'

output:
file('run_quality_report.html')

script:
"""
Rscript -e 'rmarkdown::render(input = "$rscript", output_dir ="${params.outdir}", params = list(directory = "${params.outdir}"))'
"""

Next I tried `output_dir ="$baseDir/${params.outdir}"`. And the output ended up at this path: `results/run_quality_report.html`. But still there was error `Missing output file(s)`

Does the output need to be at the `.../work/82/results/...html` path with a symbolic link from `results/...html`? Why is nextflow returning Missing output file error? How may I remedy this?

My nextflow.config parameters looks like this

params {
outdir = "../results"
}

The rscript header and some snippets looks like this

output:
html_document:
toc: true
toc_float: true
# runtime: shiny
params:
directory:
value: x
...

```{r , echo=FALSE, warning=FALSE}
inputfile <-read_tsv(file.path(params$directory, "filename.tsv"), show_col_types = FALSE)

Answer: https://www.biostars.org/p/9473324/#9473334

You just have to trick rmarkdown to spill the file to the current dir (the nextflow work dir) and Nextflow will then publish it to the outdir.
Question: I am using the module Biopython module NCBIWWW to blast some sequences online. I would like to blast my sequences against different databases available, however I cannot find a comprehensive list of them.

Here is an eample of simple query to the Nucleotide collection database using "blastn" algorithm.

from Bio.Blast import NCBIWWW

result_handle = NCBIWWW.qblast("blastn", "nt", some_sequence)

As you can see, the database Nucleotide collection is specified as "nt". With what shall I substitute "nt" in case I want to query the Human GRCh37/hg19 database for example? And if I want to query other species/builds? Is there any comprehensive list available where I can find the short names for all the databases available at http://blast.ncbi.nlm.nih.gov ?

Thanks!

Answer: This is [the list of preformatted BLAST databases available at BLAST FTP site][1].

I parsed the list and tested these databases with a script. All of them work except 3.

**Protein databases**

```
env_nr
nr
pataa
pdbaa
refseq_protein
swissprot
```

**Non-RefSeq nucleotide databases**

```
env_nt
est
est_human
est_mouse
est_others
gss
htgs
nt
patnt
pdbnt
sts
vector
wgs
```

**RefSeq nucleotide databases**

```
human_genomic #obsolete, renamed to refseq_genomic_human
human_genomic_transcript #obsolete
mouse_genomic_transcript #obsolete
other_genomic
refseq_genomic
refseq_rna
refseqgene #obsolete
```

To your second question about how to search against an older human genome build (Human GRCh37/hg19), I am afraid you cannot do it with remote BLAST. You need to download [the older data][2] and create a BLAST database and run BLAST locally.

[1]: http://www.ncbi.nlm.nih.gov/books/NBK62345/#_blast_ftp_site_The_blastdb_subdirectory_
[2]: ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/
Question: I am trying to run qualimap on multiple bam file via command.

qualimap multi-bamqc -d bam_file.txt -gff target.interval_list

But it gives an error as follows:

qualimap multi-bamqc -d bam_file.txt -gff target.interval_list
Java memory size is set to 120000M
Launching application...

OpenJDK 64-Bit Server VM warning: ignoring option MaxPermSize=1024m; support was removed in 8.0
QualiMap v.2.2.1
Built on 2016-10-03 18:14

Selected tool: multi-bamqc
Running multi-sample BAM QC

Checking input paths
Failed to run multi-bamqc
java.lang.RuntimeException: The raw data doesn't exist for sample: 001X1.sorted.deduped.bqsr.bam
Folder path:/media/molmed/Data/PJRCT1/bam
Please check raw data directory is present.

bam_file.txt conatins

`head bam_file.txt`

001X1.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
001X2.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
001M1.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
001M2.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
003X3.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
003M1.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam

I also tried by giving the complete path to the file as `/media/molmed/Data/PJRCT1/bam/001X1.sorted.deduped.bqsr.bam`

and also by adding `/` to the folder ends as `/media/molmed/Data/PJRCT1/bam/`. But all gives the same error

Answer: Not sure if it is related with copy/paste, but the file `target.interval_list` does not have any file extension.

Is this file a GFF/GTF or BED file?

Regarding the file `bam_file.txt`, the first column, I think it is the sample name that you want to give, not necessarily the same as the file name. In the second column you need to specify the bam file (this is my understanding from reading the docs - I never run this specific command), something like this:

sample_1 /media/molmed/Data/PJRCT1/bam/001X1.sorted.deduped.bqsr.bam

If you do not have `BAM QC result` you might need to specify the `-r` option:

> Raw BAM files are provided as input. If this option is activated BAM QC process first will be run for each sample, then multi-sample analysis will be performed.

[Source](http://qualimap.conesalab.org/doc_html/command_line.html#multi-sample-bam-qc) of this information.

I hope this helps,

António
Question: I am using trimmomatic for quality trimming of Paired reads. however, i am unable to get equal numbrer of reads in paired read files. Can anyone please help me out to resolve this issue. the command I am using is as follows:

java -jar trimmomatic-0.36.jar -phred33 PE /path/to/dir/SRR9590_1.fastq /path/to/dir/SRR9590_2.fastq SRR9590_1P.fq SRR9590_1U.fq SRR9590_2P.fq SRR9590_2U.fq ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

Answer: Hello,

first let's check if there are the same number of line in both fastq files. So what's the output of `wc -l /path/to/dir/SRR9590_1.fastq` and `wc -l /path/to/dir/SRR9590_2.fastq`?

If it's the same number we can go on. In your output above you have this message:

java.io.FileNotFoundException: /data/sata4/memona/Trimmomatic-0.36/TruSeq3-PE.fa

Check if this file realy doesn't exist. If there is a similar named file adopt your command line to this. I have found some threads where the authors guess that the error `thread "main" java.lang.RuntimeException: Sequence and quality length don't match:` arise due to the missing adaptor file.

fin swimmer
Question: Hi,

I have a confusion regarding the filtering of VCF file. Consider the line below:

chr1	14464	rs546169444	A	T	112.33	.	AC=1;AF=0.125;AN=8;BaseQRankSum=0.548;ClippingRankSum=0.000;DB;DP=123;ExcessHet=3.0103;FS=1.574;MLEAC=1;MLEAF=0.125;MQ=28.09;MQRankSum=0.745;QD=2.88;ReadPosRankSum=2.036;SOR=0.344	GT:AD:DP:GQ:PL	0/1:31,8:39:99:144,0,845	0/0:27,0:27:81:0,81,799	0/0:30,0:30:90:0,90,847	0/0:27,0:27:81:0,81,799

There are two DP, one in INFO field DP=123 while the other is in FORMAT field. As it has four samples so each sample has its own DP in FORMAT field: 39, 81, 30, and 27.

My confusions are:

1. The additive of DP of the FORMAT field is not 123, why is the difference???,

2. If we filter vcf file on the basis of DP, which DP would the filtering tool used, INFO's DP or FORMAT's DP?

Thanks in advance,

Waqas.

Answer: https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_annotator_Coverage.php

> "At the sample level (**FORMAT**), the DP value is the count of reads that passed the caller's internal quality control metrics (such as MAPQ > 17, for example)."

.

> "At the site level (INFO), the DP value is the **unfiltered** depth over all samples."

.

>> If we filter vcf file on the basis of DP, which DP would the filtering tool used, INFO's DP or FORMAT's DP?

depends of our needs...
Question: Dear All,

I have a data matrix that has 17 samples and over 800 genes that belong to ten different gene families. I want to show these gene families in heatmap by marking them.

I performed heatmap but I do not know how to show gene families in the heatmap graph.

Anyone knows how to that?

Answer: Take a look at the heatmap.3 and ComplexHeatmap packages to mark your genes.