cannin/biostars_qa huggingface dataset:
Question: I am newbie to linux and NGS. Can anyone help me out how to get information about splice junctions using HISAT2?

The command I'm using is giving information in single file about alignment in SAM format. The command is as follows:

./hisat2 -p 64 --max-intronlen 10000 -x /data/memona/hisat2-2.1.0/hisat_index -1 /data/memona/SRR959590_A_1P.fq -2 /data/memona/SRR959590_A_2P.fq -S /data/memona/results/hisat_align.sam &

Answer: Hi blooming.daisy333,

It is explained in the [manual][1]. You have to use the option --novel-splicesite-outfile.
**Update:**
Be careful there is an error how left splice sites are reported as I mentioned [here][2]
The splice sites are correctly reported but the format they use is not necessarily straightforward.

[1]: https://ccb.jhu.edu/software/hisat2/manual.shtml
[2]: https://github.com/infphilo/hisat2/issues/164
Question: hello all,

i am having a hard time in using bcftools set set some variants in my vcf to missing.

as i understand missing2ref is no longer used and replaced with setGT. having said that, it is not clear how to use it.

i need to set any variant with genotype quality (GQ) < 20 to missing. this will require that i use the plugin +fill-tags at the first stage, but then how do i use the setGT plugin? how is the filter command written?

i am so lost here. thanks in advance

Answer: bcftools +setGT in.vcf.gz -- -t q -n . -i 'FMT/GQ<90'
Question: I am trying to plot coverage for CDS region using deeptools.

from the manual there is a flag ``` --region ```, but it either takes a chromosome or region of chromosome

> –region chr10 or –region chr10:456700:891000

Is there is a way to pass bed file to it like in ``` –blackListFileName, -bl ``` flag?

or I need to use something like:

``` samtools view -hL cds.bed alignment.bam > alignment_in_cds.sam ```

and then use plotCoverage normally?

Answer: `plotCoverage` itself can't be made to use a BED file like you want. If you just want something like "% of signal in the CDS" then you can give a GTF file and your BAM file(s) to `plotEnrichment`, which will produce a histogram of "% alignments/feature" and (optionally) a table of the values.

If you really do want only the coverage distribution over the CDS then there's no way to do that directly with deepTools. One can prefilter the BAM file as you showed, but note that this will still create off-target bases that are flanking each CDS, so even then the results will only be approximate.
Question: Dear All,

I would like to know about few docking softwares, freely available ones.

Thanks

Answer: Autodock is one of the popular softwares available in open-source for docking. Popular commercial ones includes Tripos SYBYL, Discovery studio and schrodinger.

You can find the list of software in [wikipage][1]

[1]: http://en.wikipedia.org/wiki/Docking_%28molecular%29#Mechanics_of_docking
Question: Hi, I am trying to convert a fastq file (containing multiple sequences) to a fasta file using fastx_toolkit (Hannon lab). However I get this error message:

fastq_to_fasta: Error: invalid quality score data on line 68 (quality_tok = "ATATGCGTGCCATTG...etc)

I notice on another thread that you can put -Q33 to tell it you are using Illumina quality scores. Does anyone know if there is an equivalent flag to tell it I am using nanopore data?

If not can anyone recommend another way to convert these files?

Thank you!

Answer: You can easily do that conversion using biopython:

from Bio import SeqIO
count = SeqIO.convert("input.fastq", "fastq", "output.fasta", "fasta")
print("Converted %i records" % count)
Question: I just found the FindSubCluster tool within Seurat, and am super excited to use it. Just not sure exactly how!

The usage is here:

```
FindSubCluster(
object,
cluster,
graph.name,
subcluster.name = "sub.cluster",
resolution = 0.5,
algorithm = 1
)
```

I receive an error when using it.

> test <- FindSubCluster(scfp, 7, scfp, subcluster.name = "test", resolution = 3, algorithm = 1)
Error in match(x, table, nomatch = 0L) :
'match' requires vector arguments
> test <- FindSubCluster(scfp, 7, scfjjd, subcluster.name = "test", resolution = 3, algorithm = 1)
Error in `[[.Seurat`(object, graph.name) : object 'scfjjd' not found
> test <- FindSubCluster(scfp, 7, RNA, subcluster.name = "test", resolution = 3, algorithm = 1)
Error in `[[.Seurat`(object, graph.name) : object 'RNA' not found

My object name is scfp. I have done FindClusters already. The cluster I wish to subcluster is 7.

### I think my trouble is with the **graph.name** argument.

I tried this:

> test <- FindSubCluster(scfp, "7", "scfp", subcluster.name = "test", resolution = 3, algorithm = 1)
Error: Cannot find 'scfp' in this Seurat object

Still no success yet...

Any help would really be appreciated!

Answer: Okay so I got it... I think.
I explored the Seurat object a litle bit more and found that the cluster assignments *were* saved.
I was able to visualize using the `group.by` argument so here is the process:

scfp <- FindNeighbors(scfp, graph.name = "test", dims = 1:100)
scfp <- FindClusters(scfp, graph.name = "test", resolution = 2, algorithm = 1, verbose = TRUE)
scfp <- RunTSNE(scfp, dims = 1:100)
DimPlot(scfp, reduction = "tsne", label = TRUE, label.size = 6 )

![enter image description here][1]
And then after doing this:

scfp <- FindSubCluster(scfp, "6", "test", subcluster.name = "blood", resolution = 2, algorithm = 1)
DimPlot(scfp, reduction = "tsne", group.by = "blood", label = TRUE, label.size = 6)

This is the resulting tSNE:

![enter image description here][2]

[1]: https://i.ibb.co/W3NXrNp/before.png
[2]: https://i.ibb.co/S683P2W/after.png

Of course going to play around to optimize cluster numbers...

Thank you @fracarb8. You got me going!
Question: Hi!

I have a genelist from a KEGG pathway say TGFb signalling pathway, that has 80 genes in it (set A). and I have second list of genes with me which accounts for 45 genes (set B), and the overlap between the two sets is of 15 genes.

How can I get a p-value associated with this overlap?

Has anyone done this before?

I think, phyper of R can be applied to this .......

Thank you

Answer: Please check these questions:

- [Statistical measure of biological significance for overlapping genomic regions][1]
- [How to compute gene enrichment P-value? (significantly more than expected)][2]

[1]: /p/8480/
[2]: /p/66729/
Question: I would like to (try to) publish a software I have been developing; the software is for analyzing and integrating data from ENCODE, and of course will be open source.

I know about Bioinformatics APPLICATIONS NOTES as a possible venue, but I would like to evaluate alternatives.

Can anyone point me to some valid journal accepting this kind of publications?

Answer: Here are few journals to consider,

- BMC Bioinformatics
- Briefings in Bioinformatics
- Molecular systems biology
- PLOS Computational Biology
- GigaScience
- Nucleic Acids Research
- Genome Biology (Software)
Question: What is the difference between GENCODE and Ensembl annotation? That's actually the first question in the [GENCODE FAQs][1]:

> The GENCODE annotation is made by merging the Havana manual gene
> annotation and the Ensembl automated gene annotation. ...In practical
> terms, the GENCODE annotation is identical to the Ensembl annotation.

I am looking at the mouse data for GENCODE M15 compared to Ensembl 90, which should be comparable according to both source. Total number of transcripts is 131,100 vs 131,195, so that difference is negligible. However, some subsets are very different. The number of protein-coding genes is 21,950 vs 22,598, which is a little more noticeable. Long non-coding RNA genes is 11,975 vs 8,980, so more than 30% drop. Thus, it seems like annotation is not really identical. Are those differences real or are they just counting the gene biotypes differently?

[1]: https://www.gencodegenes.org/faq.html

Answer: The Gencode statistics webpage refers to the annotation on the reference chromosomes only, which has 131,100 transcripts.
http://www.gencodegenes.org/mouse_stats/current.html

The Ensembl statistics include all primary assembly regions, which explains the higher number of transcripts (131,195).
http://www.ensembl.org/Mus_musculus/Info/Annotation

Furthermore, the grouping of the gene biotypes for the statistics webpages differ between Gencode and Ensembl:

Gencode includes only the 21,950 genes with a "protein_coding" biotype under the "Protein-coding genes" category on the webpage. Ensembl reports in the "Coding genes" category all genes that contain an ORF, which adds the IG/TR genes for example.

In the case of the long non-coding RNA genes, the difference is due to the inclusion of the TEC genes by Gencode (~3000).
Question: Please help me find the number of mapped reads from a bam file. I was going through forums and tutorials. Looking at `samtools flagstat` resulted the following:

My total read count is **413,033**. So, I understand 453,800 to be number of alignments. Now, is it that 391,696 reads mapped? According to [**this tutorial**][1], *Yes*.

1\. `samtools flagstat file.sorted.bam`

```
453800 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
391696 + 0 mapped (86.31%:-nan%) 0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
```

Found another explanation at [this page][2]. It says 391,696 is the number of locations that reads mapped. A number same as flagstat is giving us above and the number of mapped reads are 413,033 same as my original fasta sequence number. Check commands below.

2\. Number of Mapped locations:

```
samtools view -F 0x04 -c file.sorted.bam
391696
```

3\. Number of Mapped reads:

```
samtools view -F 0x40 file.sorted.bam | cut -f1 | sort | uniq | wc -l
413033
```

Could I please know, which one, should I follow to find the number of reads that mapped the reference from my original fasta file? Suppose the reads are to be single end. I would thankful to your explanations.

[1]: http://left.subtree.org/2012/04/13/counting-the-number-of-reads-in-a-bam-file/
[2]: https://gist.github.com/davfre/8596159

Answer: 1. 391696 is the number of mapped alignments in the BAM file. That includes multimappers, so it won't necessarily be the number of mapped reads.
2. This is identical to above and does the exact same thing.
3. This is the number of entries excluding those marked as read #1 in a pair. This is not what you want.

You have a single-end data set. What you want is the number of mapped alignment, excluding those marked as secondary or supplementary (i.e., multimappers or chimeric entries). Consequently, `samtools view -F 0x904 -c file.sorted.bam` should work for you. This won't work for paired-end datasets (you'd have to `samtools view -F 0x4 foo.sorted.bam | cut -f 1 | sort | uniq | wc -l` them).
Question: Hi everyone,

I have a microarray set of 26933 genes and 6 samples (3 controls and 3 mutants). I calculated differential expressed genes using limma's topTable() function. I got a result of this format:

> head(res.new)

logFC AveExpr t P.Value adj.P.Val B
17291881 -2.157936 9.021125 -41.49770 3.218440e-09 3.829474e-05 10.82552
17357688 -2.364349 7.160051 -41.15755 3.397181e-09 3.829474e-05 10.79753
17224251 1.820522 10.665778 37.01446 6.816679e-09 3.829474e-05 10.41426
17239845 -2.025497 10.054660 -35.87667 8.366167e-09 3.829474e-05 10.29350
17304523 -1.830077 10.117184 -35.36567 9.191590e-09 3.829474e-05 10.23680
17229782 2.716749 7.911040 35.24054 9.407763e-09 3.829474e-05 10.22267

The logFC values are log2(FC). How do I convert the values (positive and negative) to normal fold changes?

I tried this:

#take 2^(absolute(log2FoldChange))

> res.new$FC<-2^(abs(res.new$logFC))

#change the sign of FC according to log2FC

> res.new$FC<-ifelse(res.new$logFC<0,res.new$FC*(-1),res.new$FC*1)

But I am not sure if I am correct. Please let me know.

Thanks!

Answer: res.new$FC <- 2^res.new$logFC

Negative log2 fold-changes are just less than 1.
Question: Hello,

I have a deep sequencing experiment to analyze and I am hesitant about the variant caller algorithm/program to use as I have some doubts concerning scalability

For a small experiment, not a deep-seq one, we usually rely on the gatk recommendations and best practicies guides, by combining a couple of tasks such as deduplication, this makes the analysis time more or less acceptable. For a deep sequencing experiment, there is no rational about removing duplicates, which makes the variant call part very long and relatively not scalable.

I didnt try the gatk variant caller yet on deep-seq data but I guess it will take a lot of time, I was wondering what is the best option to do the variant call on a deep seq data in term of scalability, if anyone has tried that in the past, I would be grateful if he can give some hints on the best way of doing this.

Thanks

Answer: For 90x coverage, GATK works fine for us. You can parallelize GATK by running it per chromosome. If you want to go even further, freebayes can be safely parallelized to non-overlapping regions of chromosomes. In any case, for most callers, 90x shouldn't be too bad. Do make sure that any high-depth filters that are inherent in the defaults are either turned off or set to more sane numbers for your data.
Question: Dear all, please could you advise :

what is the best strategy to call the somatic variants in the absence of germline sample (when only the tumor sample is available) ? thank you,

-- bogdan

Answer: If your goal is to discover new somatic variants, I am afraid that might be quite challenging because the VAST majority of variants in any sample will actually be polymorphisms and not somatic variants. Use of a set of databases of "normal" variation to remove as many germline variants as possible can reduce the false positive results, but this approach will be clearly inferior to having a matched normal. If, on the other hand, the goal is to identify known cancer variants present in the samples being sequenced, having tumor only is not at all unusual and is a viable strategy.
Question: Hola!

Recently, I heard a talk about using machine learning in ChIP-Seq data analysis and binding site predictions, thus got interested. I wasn't able to find some solid reviews or examples on, what one can do and achieve with it. How powerful is it in this context and are the results reliable enough. Could you shed some light on

the aspect of using machine learning with regard to ChIP-Seq (and/or RNA-Seq),

what to expect and

the advantages/disadvantages, last but not the least,

what/how to start about it.

Few helpful links :

Application of machine learning methods to histone methylation ChIP Seq data reveals H4R3me2 globally represses gene expression

Integration of ChIP-seq and machine learning reveals enhancers and a predictive regulatory sequence vocabulary in melanocytes.

Thanks

Answer: Here are some pointers on transcription factor binding prediction and Chip-seq methods. Most of them are not strictly machine learning papers, but one or the other method of ML has been employed to discover "new biology".

Predicting binding from sequence

Weirauch, M.T. et al., 2013. Evaluation of methods for modeling transcription factor sequence specificity. Nature Biotechnology.

Hannenhalli, S., 2008. Eukaryotic transcription factor binding sites-modeling and integrative search methods. Bioinformatics, 24(11), pp.1325–1331.

Zambelli, F., Pesole, G. & Pavesi, G., 2009. Pscan: finding over-represented transcription factor binding site motifs in sequences from co-regulated or co-expressed genes. Nucleic Acids Research, 37(Web Server), pp.W247–W252.

ChIP-seq based

Gerstein, M.B. et al., 2012. Architecture of the human regulatory network derived from ENCODE data. Nature, 489(7414), pp.91–100.

Wang, J. et al., 2012. Sequence features and chromatin structure around the genomic regions bound by 119 human transcription factors. Genome Research, 22(9), pp.1798–1812.

Whitfield, T.W. et al., 2012. Functional analysis of transcription factor binding sites in human promoters. Genome Biology, 13(9), p.R50.

Abstracting tissue specificity by DNAse hypersensitivity

Arvey, A. et al., 2012. Sequence and chromatin determinants of cell-type-specific transcription factor binding. Genome Research, 22(9), pp.1723–1734.

Natarajan, A. et al., 2012. Predicting cell-type-specific gene expression from regions of open chromatin. Genome Research, 22(9), pp.1711–1722.

Pique-Regi, R. et al., 2011. Accurate inference of transcription factor binding from DNA sequence and chromatin accessibility data. Genome Research, 21(3), pp.447–455.
Question: The SAM specification states the following:

> BGZF is block compression implemented on top of the standard gzip file format....a compliant gunzip utility can decompress a BGZF compressed file

I would like be able to read a BAM file using gzip, as the specification says I should be able to. However, when I open a BAM file using the inbuilt python library gzip I am able to read the header normally but the reads still look compressed in some manner, like so:

```
AMCBCZ0XDZ100SMiASi??'"?Icd(ZHS2000-943_70:2:2301:14923:121232@("("("("("("("("("("("("("("("("(""%%##%'&##'((&''($'((($'($&%&'%%&&'(("'#%#&'&&$&(&(( !%!ASi?'!" !!!!""BCZ0XDZ100SMiIAd?sHS2000-905_68:3:1102:3166:478590?B ?"?"?"?"?"?"?"?"?"("("("("("("("!"%%##%''''#''!'(&'%'((%$&&(('('"'!'&(&!'&"# %&#&&$"!! "BCZ0XDZ51^1$30^CCCT$18SMiASi?'"I?d'HS2000-943_70:2:1302:17909:144833"("("?"?"?"?"?"?"?"("(!"!"?$!(!!"!"!"%%# #'#%!#& #%$###%%"!'(#'"!H0CH1CH2CBCZ0)'"IYd????????HS2000-943_70:2:1302:17909:144833@($""(!?"?"$?"?"?"??""("("("("("("("!
```

I can see the making of a real read in there, (i.e the read name: `HS2000-943_70:2:2301:14923:121232`) but most of it is illegible. Perhaps I should be reading in entire blocks and then decompressing further? I am not sure. I have also tried reading the file using the bgzf biopython module and receive the same results.

I am aware of [the BioStars question here][1]. However, I feel given the explicit declaration in the SAM file specification that the files are readable by gzip, the answers (to use pysam) was inadequate.

Thanks for reading!

[1]: https://www.biostars.org/p/78926/

Answer: Well the reads are still in binary, rather than in a text format like in a SAM file (i.e., decompressing a BAM file won't give you a SAM file). The format for a BAM file is mentioned in the SAM spec., so you'll need to write an interpreter to parse the binary (it'd be much easier to just use pysam).
Question: I downloaded the gff3 annotation file for [Zebrafish][1] from NCBI and ran [gff3validator][2] from genometools as such

gt gff3validator GCF_000002035.5_GRCz10_genomic.gff

and I get the following error:

gt gff3validator: error: CDS feature on line 621884 in file "GCF_000002035.5_GRCz10_genomic.gff" has the wrong phase 1 (should be 0)

I have not modified the file in any way. What does this error mean? Should I be concerned that it could lead to potential issues? How can this be fixed? Thanks.

[1]: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_other/Danio_rerio/latest_assembly_versions/GCF_000002035.5_GRCz10
[2]: http://www.genometools.org/tools/gt_gff3validator.html

Answer: I have become interested in tracking down this and ran what the OP did. Basically I am trying to determine if this is a GFF file format error or a validation error.

I do get the same error as the OP:

gt gff3validator: error: CDS feature on line 621884 in file "GCF_000002035.5_GRCz10_genomic.gff" has the wrong phase 1 (should be 0

The obtain the line:

cat GCF_000002035.5_GRCz10_genomic.gff | awk ' NR==621884 { print $0 } '

This produces the offending line:

NC_007121.6	BestRefSeq	CDS	21750127	21750129	.	+	1	ID=cds20891;Parent=rna31380;Dbxref=GeneID:553997,Genbank:NP_001019280.1,ZFIN:ZDB-GENE-050609-5;Name=NP_001019280.1;Note=The RefSeq protein has 9 substitutions compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=CDS;gene=pcdh1g9;product=protocadherin 1 gamma 9;protein_id=NP_001019280.1

Looking at the genbank format at https://www.ncbi.nlm.nih.gov/nuccore/NC_007121.6 it shows:

join(21744267..21746679,21746681,21750127..21750129,
21750131..21750164,21855170..21855228,21858956..21859074,
21860119..21860155,21860973..21860979)
/gene="pcdh1g9"
/gene_synonym="DrPcdh1g8"
/inference="similar to AA sequence (same
species):RefSeq:NP_001019280.1"
/exception="annotated by transcript or proteomic data"
/note="The RefSeq protein has 9 substitutions compared to
this genomic sequence; Derived by automated computational
analysis using gene prediction method: BestRefSeq."
/codon_start=1
/product="protocadherin 1 gamma 9"
/protein_id="NP_001019280.1"
/db_xref="GI:66773380"
/db_xref="GeneID:553997"
/db_xref="ZFIN:ZDB-GENE-050609-5"

This shows that it is the third CDS that raises the error. Add up the lengths of the previous CDS sizes and see how far are we from the multiple of 3. That would be the phase.

>>> size = 21746679 - 21744267 + 1 + 1
>>> divmod(size, 3)
(804, 2)

The remainder is 2, this means that the next codon starts one base in. So phase should be 1.

Basically telling us that the GFF is correct and the validator is incorrect.
Question: I am using the enricher function of clusterprofiler for GO enrichment of my non-model organism.

First I have my DE genes (n=638, converted as character):

gene <- read.csv("list.csv", header = F,sep=",")
gene <- as.character(gene[,1])
head(gene)
[1] "BANY.1.2.t20473" "BANY.1.2.t12787" "BANY.1.2.t10473" "BANY.1.2.t10472"
[5] "BANY.1.2.t08098" "BANY.1.2.t04432"

Then I have the term2gene:

term2gene <- read.csv("GO2gene.csv",header=F,sep=",")
head(term2gene)
V1 V2
1 GO:0002119 BANY.1.2.t15918
2 GO:0003916 BANY.1.2.t01379
3 GO:0005344 BANY.1.2.t10849
4 GO:0005344 BANY.1.2.t10852
5 GO:0005344 BANY.1.2.t10851
6 GO:0005344 BANY.1.2.t11738
and also term2name:

term2name <- read.csv("GO2name.csv",header=F,sep=",")
head(term2name)
V1 V2
1 GO:0000001 mitochondrion inheritance
2 GO:0000002 mitochondrial genome maintenance
3 GO:0000003 reproduction
4 GO:0000006 high-affinity zinc transmembrane transporter activity
5 GO:0000007 low-affinity zinc ion transmembrane transporter activity
6 GO:0000009 alpha-1,6-mannosyltransferase activity

Then run it:

x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.05, pAdjustMethod = "BH",qvalueCutoff = 0.1)

And I got this...

--> No gene can be mapped....
--> Expected input gene ID: BANY.1.2.t09827, BANY.1.2.t01112, BANY.1.2.t16389, BANY.1.2.t23551, BANY.1.2.t13938, BANY.1.2.t11911
--> return NULL...

Everything seems fine so I really hope to know what the problem is, thank you.

Answer: Hello, can you run `intersect(term2gene$V2, gene)` ? Maybe no gene hit by any pathway.