cannin/biostars_qa huggingface dataset:
Question: Hi,

I'm new to the NGS data analysis. i am studying from the very basics of NGS how the data are analysed.

For aligning the Raw reads there many tools each one have different algorithm for
example
BWA works on BWT
Novoalign works on Needleman wunch
Mosaik works on Smith- Water man

what about the variant caller
Most of the people are using GATK, SAMtools mpileup, and freebayes for the variant calling i have gone through their site but i cant find what algorithm are the using for detecting the variant where in freebayes they told that they have been using the bayesian. i want to know whether all the variant caller works on the bayesian algoritm or they are having different way of finding the variant using different algorithm

Answer: I found these are the major difference between GATK and Samtools

1. Preprocess of alignment" GATK drops reads with low mapping quality , but samtools uses all reads by default.
2. SNP genotype likelihood model: GATK assumes sequencing errors are independent while samtools believes the second error comes at a higher chance.Indel genotype likelihood model. This is one of the major differences between samtools and GATK. Samtools model is derived from BAQ. GATK's model was derived from Dindel's model.
3. Samtools uses hand-tuned filters, while GATK learns filters from data. Of course GATK's approach is more convenient and powerful at least for human variant calling where you have enough data to train the model and you do not need non-polymorphic sites
4. Other GATK specific features: haplotype caller, indel realignment among many others.
5. Other samtools specific features: genotype-free analysis, physical phasing and so on.
Question: Does anyone know where I can get the chrom sizes info for the ensembl reference genome. I'm attempting to do bedToBigBed but need the chrom sizes file and can't seem to find the ensembl version.

EDIT: sorry, should have been more clear on that. I used the [igenomes package](https://support.illumina.com/sequencing/sequencing_software/igenome.ilmn) for Human (Grch37) - To run the Tuexedo pipeline. I got my GTF file of all the merged transcripts out, converted that to BED for use in an web based genome browser I'm working on. The required binary format is bigbed, so I'm using the bedToBigBed utility from UCSC which requires a chrom info file (start points of each chromosome + patches in the reference Fasta sequence) in tab delimited text format.

UPDATE: So, here's how I solved my problem in case anyone else comes across this.... I basically got rid of all the patches in the GTF file. That can be done with awk (I'm sure there are more efficient methods using regex https://gist.github.com/AndrewSkelton/10829401). The reference genome patches, as far as I can tell, are not utilised by any of the Tuexedo packages and just get in the way when doing file conversions and such in downstream. The solutions below using samtools and using looking in the SAM headers give the chromosonal lengths, but not the patch lengths. Thanks for all the solutions!

Answer: So, here's how I solved my problem in case anyone else comes across this.... I basically got rid of all the patches in the GTF file. That can be done with awk (I'm sure there are more efficient methods using regex https://gist.github.com/AndrewSkelton/10829401). The reference genome patches, as far as I can tell, are not utilised by any of the Tuexedo packages and just get in the way when doing file conversions and such in downstream. The solutions below using samtools and using looking in the SAM headers give the chromosonal lengths, but not the patch lengths. Thanks for all the solutions!
Question: QDQSDCCSDQSDQQCQSCSDQSSQDCQSQSCQDQSDQSDQSD

I followed this link from stackoverflow to place a marker on google maps. Currently i am doing like this to open google maps using iOS app swift.

UIApplication.shared.open(URL(string:"comgooglemaps://?center=\(self.location.coordinate.latitude),\(self.location.coordinate.longitude)&zoom=14&q=\(self.location.coordinate.latitude),\(self.location.coordinate.longitude)")!, options: [:], completionHandler: nil)
its opening map and placed a marker which is fine, but my problem is its showing lat/lng on marker, and i want to add description (which i already have to display) on marker instead of lat/lng.

Answer: public function teacherRecords(){
$teachers= DB::table('teachers')
->leftjoin('religions', 'teachers.religion_id', '=', 'religions.id')
->leftjoin('areas', 'teachers.area_id', '=', 'areas.id')
->select('teachers.*')
->where('teachers.id',$this->id)
->first();
}
Question: I want to break my bed file into the chromosome-specific files containing the data of respective chromosome so I used the following command

for chr in {1..22}; do awk -v chr=$chr '$1 == chr' data.txt > data_chr$chr.txt; done

The issue is that the files created do not contain any information.

My bed file has six columns and is having the following format

chr1	10094	10144	563_280_1408	1	+
chr1	10095	10144	421_286_567	0	+
chr1	10165	10215	382_570_1514	9	+
chr1	10170	10196	664_789_434	1	+
chr1	10171	10221	439_96_1067	10	+
chr1	10174	10223	378_525_47	5	+
chr1	10174	10224	94_430_1386	6	+
chr1	10217	10279	4_218_1107	5	+
chr1	10224	10267	138_1957_415	59	-
chr1	10242	10272	638_1023_1664	4	+

Any suggestions will be highly appreciated

Answer: You need to prefix the `chr` variable with the string literal `chr`, as right now you are asking `awk` to match the first field on numbers (1 through 22) and not chromosome names (chr1 through chr22).

$ for chr in {1..22}; do awk -v chr="chr${chr}" ...

But BEDOPS offers a much faster and simpler way of splitting BED files with [bedextract][1], which uses a binary search approach on sorted BED files to jump to the start of a chromosome:

$ fn=intervals.bed
$ for chr in `bedextract --list-chr ${fn}`; do bedextract ${chr} ${fn} > ${fn}.${chr}.bed; done

With `awk` or `grep`, you need to re-read through the file linearly on every chromosome, which wastes time, especially on very large, whole-genome scale BED files. Use BEDOPS if your time is valuable to you.

[1]: http://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedextract.html
Question: I have found the `bigPint` package in bioconductor recently and found it very helpful to plot the data for QC reasons.
The vignettes are very good made and the package is very neat, but unfortunately they are all working with `edgeR`.
Even though it says they can also run with a `DESeq` object, there is no example of such an analysis.

I was wondering if the authors of this package can also add an example with `DESeq2` workflow. Or at least mention what columns of the `deseq` object are needed in creating the plots.

thanks

Answer: Hello [Assa Yeroslaviz][1]. Thank you for this suggestion. We added reproducible example analyses and more detailed explanations of columns using `DESeq2`. You can see more information [here][2].

[1]: https://www.biostars.org/u/2295/
[2]: https://github.com/lindsayrutter/bigPint/issues/5
Question: Hi,

I am relatively new to bioinformatics/R/scripting. Currently I am trying to improve my coding style so that my scripts are more accessible, both for me and other readers.

Does anyone know of any resources for model scripts I could use as inspiration?

Thanks.

Answer: **This will be good starter for you,**

https://cran.r-project.org/doc/contrib/Paradis-rdebuts_en.pdf (I have started with this PDF :-))

http://www.tutorialspoint.com/r/

**Advanced usage of available pipelines:**

For some high throughput analysis R scripts/[packages][1] look into [Bioconductor][2] website.

**More discussions,**

https://www.biostars.org/p/539/

[1]: https://bioconductor.org/packages/release/BiocViews.html#___Software
[2]: https://bioconductor.org/
Question: Hi everyone,

I am trying to parse a table file I obtained with hmmscan by comparing a set of fastas against the whole Pfam database of protein domains. I am trying to use HmmerIO in Python to parse the ouput, but it seems to me that, for each fasta, it only keeps the domain with the highest evalue, while I want to keep all of those with a significant evalue.

Example:

Input file:

# --- full sequence ---- --- best 1 domain ---- --- domain number estimation ----
# target name accession query name accession E-value score bias E-value score bias exp reg clu ov env dom rep inc description of target
#------------------- ---------- -------------------- ---------- --------- ------ ----- --------- ------ ----- --- --- --- --- --- --- --- --- ---------------------
GST_N PF02798.17 gi|528981796|ref|NC_021285.1|_5567 - 4.1e-16 58.9 0.0 6.9e-16 58.2 0.0 1.4 1 0 0 1 1 1 1 Glutathione S-transferase, N-terminal domain
GST_N_2 PF13409.3 gi|528981796|ref|NC_021285.1|_5567 - 3.3e-14 52.8 0.0 5.8e-14 52.0 0.0 1.4 1 0 0 1 1 1 1 Glutathione S-transferase, N-terminal domain
GST_N_3 PF13417.3 gi|528981796|ref|NC_021285.1|_5567 - 3.4e-14 52.8 0.0 5.4e-14 52.2 0.0 1.3 1 0 0 1 1 1 1 Glutathione S-transferase, N-terminal domain
GST_C PF00043.22 gi|528981796|ref|NC_021285.1|_5567 - 7e-12 45.3 0.0 9.6e-12 44.8 0.0 1.2 1 0 0 1 1 1 1 Glutathione S-transferase, C-terminal domain
GST_C_2 PF13410.3 gi|528981796|ref|NC_021285.1|_5567 - 6.1e-07 29.2 0.4 9.3e-07 28.6 0.4 1.3 1 0 0 1 1 1 1 Glutathione S-transferase, C-terminal domain
GST_C_3 PF14497.3 gi|528981796|ref|NC_021285.1|_5567 - 0.00011 22.2 0.0 0.00017 21.6 0.0 1.2 1 0 0 1 1 1 1 Glutathione S-transferase, C-terminal domain

Here I found several statistically significant domains for the PF02798.17 gi|528981796|ref|NC_021285.1|_5567 query (fasta sequence), but the SearchIO output I have been able to get looks like this:

for qresult in SearchIO.parse(file1, 'hmmer3-tab'):
print qresult.hits[0]

Query: gi|528981796|ref|NC_021285.1|_5567

Hit: GST_N
Glutathione S-transferase, N-terminal domain
HSPs: ---- -------- --------- ------ --------------- ---------------------
# E-value Bit score Span Query range Hit range
---- -------- --------- ------ --------------- ---------------------
0 6.9e-16 58.20 ? [None:None] [None:None]

I have no idea of how to access the other hits, even if I am quite sure there has to be a way.

Anyone has any experience with this?

Thanks a lot.

Answer: I found the answer myself, it was a dumb coding mistake...
If anyone ever needs it, the point is that every fasta sequence can have different hits and hits[0] will only show the first of the list.

The correct way to proceed is:

for qresult in SearchIO.parse(file1, 'hmmer3-tab'):
for item in qresult.hits:
print item
print item.evalue # this will print the evalue of each domain hit
Question: Hi all,

I have expression measure of genes (two sets: target and contro genes) as measured by RNAseq, which is computed in tabular matrix.

head(A)
V1 V2 V3 V4 V5 V6
1 chr20:38788661-38812050_ENSDART00000136771 1.912 1.912 1.912 4.770 0.703
2 chr4:13971999-13980660_ENSDART00000067040 0.502 0.502 0.502 1.242 0.011

V7 V8 V9 V10 V11 V12 V13
1 0.137 6.534 4.689 7.086 2.646 4.376 0.684
2 0.007 0.684 1.125 2.819 2.392 1.406 0.364
head(B)
V1 V2 V3 V4 V5 V6
1 chr19:24974067-24977991_ENSDART00000052421 0.005 0.005 0.005 0.025 0.139
2 chr8:26148102-26150871_ENSDART00000049793 2.312 2.312 2.312 6.153 1.109
3 chr8:47908733-47913950_ENSDART00000025620 2.446 2.446 2.446 9.660 5.411
V7 V8 V9 V10 V11 V12 V13
1 0.807 0.010 0.232 5.963 0.022 0.149 3.162
2 0.291 4.556 3.214 16.486 2.158 5.599 0.939
3 7.745 11.614 13.150 51.680 8.557 18.356 36.378

PLOT:

par (mfrow=c(1,2))
plot (log2(A$V2/A$V5), log2(A$V8/A$V11), xlab="log2(FC)", ylab="log2(FC)", xlim=c(-5,5), ylim=c(-5,5), col="red")
abline(v=0, lty=2); abline (h=0, lty=2)
plot (log2(B$V2/B$V5), log2(B$V8/B$V11), xlab="log2(FC)", ylab="log2(FC)", col="blue", xlim=c(-5,5), ylim=c(-5,5))
abline(v=0, lty=2); abline (h=0, lty=2)

How could i plot both in the same plot?

cheers

Chirag

Answer: If I were you I would try using your first plot as is, but for the second plot instead of using the plot function, use points e.g.

plot (log2(A$V2/A$V5), log2(A$V8/A$V11), xlab="log2(FC)", ylab="log2(FC)", xlim=c(-5,5), ylim=c(-5,5), col="red")
points(log2(B$V2/B$V5), log2(B$V8/B$V11), col="blue")
abline(v=0, lty=2); abline (h=0, lty=2)

You could also use `par(new=TRUE)` instead. i.e.

plot (log2(A$V2/A$V5), log2(A$V8/A$V11), xlab="log2(FC)", ylab="log2(FC)", xlim=c(-5,5), ylim=c(-5,5), col="red")
par(new =TRUE)
plot (log2(B$V2/B$V5), log2(B$V8/B$V11), xlab="", ylab="", col="blue", xlim=c(-5,5), ylim=c(-5,5))
abline(v=0, lty=2); abline (h=0, lty=2)

edit:

Also, if the `A$V2` etc doesn't work you may want to use `A[,2]` instead.
Question: I have some older VCF files that don't have the contig length set in the VCF header. This means that Picard and some other tools that are very strict with the VCF spec won't accept them.

The contig entries in the header should be

##contig=
##contig=

but they are

##contig=
##contig=

I know that I can manually fix this by doing the following steps.

- unzipping the file
- extracting the header
- lookup the contig lenghts in a fasta.fa.fai file
- adding the lenght to the contigs records in the header with vim
- re-header with bcftools
- bgzip and tabix the the re-headered vcf file

In my hands this works but if you make a slight VIM copy paste error you will have spend a lot of time reheadering and bgzipping a large VCF file for nothing.

Therefore I would like to have a more robust automatic solution where I just give the VCF file and the reference genome file and the header is automatically fixed and a new bgzipped VCF file is written out.

Is there a tool that can add the contig lenghts to the VCF header and write out a new bgzipped VCF file?

Answer: Use awk to insert the contig before the #CHROM line. Something like

awk '/^#CHROM/ { printf("##contig=\n##contig=\n");} {print;}' in.vcf > out.vcf

you can generate the printf line above from the ref.fai file:

awk '{printf("##contig=\\n",$1,$2);}' ref.fai
Question: I used this piece of biopython codes written by Eric Normandeau to parse several sequences from a Fasta file. The codes worked nicely. I would like to have my parsed output from this format

>DR179241 similar to UniRef100_A5HIY3 Cluster: Thaumatin-like protein; n=1; ...... ATAATCTTTAGATCAGTCATCAATCTCAACAGTATCGCTTTCAATTCTCTTTCATATTGC ATGGAAGTGTGTAAATACAATTAGGGCATTCATTGAGTTGACTTCATTTAAGCGCT......

exactly like this format:

>DR179241
ATAATCTTTAGATCAGTCATCAATCTCAACAGTATCGCTTTCAATTCTCTTTCATATTGC
ATGGAAGTGTGTAAATACAATTAGGGCATTCATTGAGTTGACTTCATTTAAGCGCT......

Could anyone kindly please light me for example how to modify this piece of biopython codes written by Eric Normandeau in order to suit my purpose?

Thank you very much and have a nice day.

Answer: How about this (untested):

#!/usr/bin/env python
import sys
from Bio import SeqIO

fasta_file = sys.argv[1] # Input fasta file
result_file = sys.argv[2] # Output fasta file

def modified(records):
for record in records:
#Clear the description
record.description=""
yield record

records = modified(SeqIO.parse(fasta_file),"fasta")
count = SeqIO.write(records, result_file, "fasta")
print "Converted %i records" % count

This used a generator function to modify the records (memory efficient - you should avoid loading everything into memory). The key point for your request was to set the description to an empty string.
Question: Hi everyone,

I am a PhD student in biochemistry, and I am learning about gene expression signature. My lab generated a 36-gene mouse signature. These genes are all highly expressed. I am interested in identifying "mouse-like" human samples from a large set of primary breast tumors.

I was wondering if someone could please give me the general guidelines on how to apply a gene signature. I can write code, but don't understand the principles (I am reading tough). Is it based on the gene names and their fold-change or just the names? I am sorry for asking such a basic question, but I am learning this aspect of bioinformatics. I read that a naive Bayes classifier is a good idea? Alternatively, ranking the samples (based on how well they express the signature) and using bootstrap resampling?

I would also greatly appreciate to be redirected to a former post or tutorial.

Thank you!

Answer: One possibility (that wouldn't even require writing much code) is to use this signature as a gene set and use GSEA on the human samples to look for samples in which that set is more highly expressed than expected. The general idea is to perform GSEA on a large number of samples, many of which you expect to not show enrichment, and then look at the resulting enrichment score (or p-value) distribution. From that, you should be able to get an idea of whether the expression of this set generally follows a normal distribution or whether there's a bimodal distribution...meaning that there's a subset of samples that you're going to be very interested in. You could alternatively use resampling there, though I think it'll be quicker and easier to just have a look at the distributions first (nothing is preventing you from doing both).

That's one fairly straight forward possibility, though there are others.
Question: Hi,

I noted that some of my genomes are names with an NZ_ and other with NC_ identifier. What is the difference between them? Specifically, are NZ_ "complete" genomes?

Answer: From [The NCBI Handbook, ch18][1]: (via [this answer][2])

**NC** - "Complete genomic molecule, usually reference assembly"

**NZ** - "Unfinished WGS"

[1]: http://www.ncbi.nlm.nih.gov/books/NBK21091/table/ch18.T.refseq_accession_numbers_and_mole/
[2]: https://www.biostars.org/p/105529/
Question: Hi,

I have recently shifted from GATK 3 to GATK4. It has changed very much. I figured most of the step but one thing is giving me problem is that how to set temporary directory option in GATK4. I have used this option `export _JAVA_OPTIONS=-Djava.io.tmpdir=./tmp ` with this option `--TMP_DIR ./tmp/`. As mentioned in [this link][1] But it is giving me this error

A USER ERROR has occurred: 'export' is not a valid command.

Can any one please tell me how to resolve this error.
[1]: https://gatkforums.broadinstitute.org/gatk/discussion/comment/46418#Comment_46418

Answer: You are probably executing `gatk` as:

gatk export _JAVA_OPTIONS=-Djava.io.tmpdir=./tmp

What you want, probably (and not tested), is:

export _JAVA_OPTIONS=-Djava.io.tmpdir=./tmp

gatk

Or just (what I would do):

gatk --TMP_DIR=./tmp
Question: I tried to create a function to extract `SNP` from `fasta` file witch content alignments sequences.
Well, I create a function that make statistics for a list of N aligned sequences.
Using `Rcpp api`, the code is:

// [[Rcpp::export]]
NumericMatrix alin_stat(CharacterVector &alin,int &N_sequence)
{
int sequence_size=alin[0].size();
NumericMatrix stat_m (5,sequence_size); // set at 0 by default

for(int i=0;i Note that:
>
> alin contains aligned sequences
>
>
> N_sequence is number of sequences

Question:

- Witch condition I have to rich to select an SNP?
- Is the alignment procedure (parameters) affects the target result?

Answer: I am not sure whether this question is related to R (as stated in the tags)- the only apparent relation is that your statistics function is coded using Rcpp. For you first question, I am assuming you are asking for methods for [SNP calling][1]. If that is correct, [this older question and related answers][2] in this site might be relevant. For your second question the answer is yes, if you change the alignment parameters you may end up with a different alignment. If you end up with a different alignment then you might end up with different genomic positions called as SNPs. **Disclaimer:** I have no direct experience working with SNP data.

**EDIT**

In addition to my original response, and if you are working in R, I would take a look at the [snpStats package][3].

[1]: https://www.biostars.org/p/944/
[2]: https://www.biostars.org/p/8237/
[3]: http://bioconductor.org/packages/release/bioc/html/snpStats.html
Question: How to make a TSS enrichment plot in ATAC-seq analysis, which shows the TSS enrichment of nucleosome-free fragments, mono-nucleosome fragments, and di-nucleosome fragments?

![enter image description here][1]

I see a tutorial (https://rockefelleruniversity.github.io/RU_ATAC_Workshop.html) about TSS enrichment plot, which uses `soGGi::plotRegion` function, I am wondering what is meaning of `score` in the plot? How should I normalize the score value based on the total reads?

![enter image description here][2]

[1]: /media/images/01ecf3cd-a185-4ae5-bcbe-96b96d9e
[2]: /media/images/386e3b84-2ef9-4c35-8c0f-8603af09

Thanks

Answer: Not sure what the Score is actually quantifying, but you can do a simple min-max normalization to force the score values to be between 0 and 1 (ie fraction of signal). This looks like a glorified wrapped of the `rtracklayer` and `GViz` packages.

Truthfully, I found using the `Deeptools` suite really easy. Customization might be a little trickier, but it was super fast and the raw data is output so you can bring that into R and plot using `ggplot` if that's more your speed.
Question: I have a list of of names of molecules, some of them rather free-form (e.g "Thiolated L-alanine"). I would like to convert these names into SMILES strings. Is there any tool that could assist with dealing with these ambiguous names?

Answer: [The Pub Chem Identifier Exchange Service][1] will certainly assist you for the more standardised names that map to CIDs (that you can download as SMILES)

You can also try Chemicalize.org (if the SMILES download is working)

However, there is (by definition) no solution to ambiguous and/or non-standard names (free form as you put it). You will have to eyeball - or better still, find the real person who has invented the free form terms you are having to deal with and get them to substitute clean standard ones such as IUPAC names (e.g. Thiolated L-alanine is even google -ve)

[1]: https://pubchem.ncbi.nlm.nih.gov/idexchange/idexchange-help.html
Question: Hello,

I'm having trouble obtaining the exact chromosome number, when using standalone BLAST. My current setting is like below:

```
/...BLAST/blast-2.2.29+/bin/blastn \
-db GPIPE/9606/105/GCF_000001405.25_top_level \
-query test.fasta \
-out test.out \
-task blastn -outfmt 6 \
-remote \
-word_size=15
```

where the database is GRCh37.

The result looks like:

sample1 gi|441431411|ref|NW_004166863.1| 100.00 41 0 0 38 78 527422 527382 1e-012 77.0

As you can see, there is gene id and the position of the reference sequence, but there is no chromosome number. I couldn't find this option in the format specifiers blast provided either. I know I could get what I want by standalone BLAT, but I need to ensure BLAST also works, because it has a windows version. I also understand if I use a XML or query-anchored format, I could locate the chromosome information in a hard way. However I'm just disappointed that I can't find it in the tabular format... Could some help me please?

Thank you!

Helene

Answer: Is the chromosome number in the full name of the gene 'gi|441431411|ref|NW_004166863.1|'?

If it is, you can customize the output of blast+, the BLAST user manual has the full options here: http://www.ncbi.nlm.nih.gov/books/NBK1763/#_CmdLineAppsManual_Appendix_C_Options_for_

In your case, you might want to try

```
/...BLAST/blast-2.2.29+/bin/blastn \
-db GPIPE/9606/105/GCF_000001405.25_top_level \
-query test.fasta \
-out test.out \
-task blastn \
-outfmt "6 qseqid salltitles pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
-remote \
-word_size=15
```

Most of this is the standard format, except that I have replaced 'sseqid' by 'salltitles', which gives you all possible titles for that gene instead of just the ID. Maybe that will print your chromosome name?

The linked above manual has all possible options, so you can trim it down more depending on what you want.

There is a longer explanation in this blog post: http://blastedbio.blogspot.com.au/2012/05/blast-tabular-missing-descriptions.html The interesting parts start after 'Update (1 April 2013)'