cannin/biostars_qa huggingface dataset:
Question: I'm comparing a control vs Treatment A and control vs.Treatment B with DESeq2 though the pvalues and padj are really high. Most padj are 0.99 or so, and most pvalues range anywhere between 0.1 to 0.9. I'm suspecting taht the way I wrote the chunk of code below doesn't compare well C vs A and C vs B, what's wrong in the following code?..

```
# DESeq1 libraries
library( "DESeq2" )
library("Biobase")

r12 = read.csv("20140807_R12.csv", header = TRUE, row.names=1)

samples <- data.frame(row.names=c("C", "A", "B"), condition=as.factor(c("C", "A", "B")))
samples$condition <- relevel(samples$condition, "C")

dds <- DESeqDataSetFromMatrix(countData = as.matrix(r12), colData=samples, design=~condition)
dds$condition <- factor(dds$condition, levels =c("C", "A", "B"))

dds <- DESeq(dds)
res <- results(dds)
head(res)

DataFrame with 6 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj

3BHSD 1.946289 -0.13137862 0.3349360 -0.3922499 0.6948736 0.9963381
5_8S_rRNA 25.628958 0.05835918 0.4797943 0.1216338 0.9031891 0.9963381
5HT2A 2.748903 -0.08777809 0.3950686 -0.2221844 0.8241703 0.9963381
5S_rRNA 1179.943635 0.39496852 0.4627465 0.8535312 0.3933648 0.9963381
7SK 1325.918947 0.39695551 0.4538413 0.8746571 0.3817605 0.9963381
A1BG 5.639654 -0.28433104 0.4541317 -0.6260982 0.5312506 0.9963381
```

Answer: There's no fix because there's nothing to fix. I'm a little surprised you got adjusted p-values that low even. There are no meaningful p-values for your experiment because computing them would be impossible. No replicates means no meaningful p-values.
Question: Hello Everyone,

I was trying freebayes on a small bam file and wanted to compare the quality of its output -vcf1- against the quality of the output of gatk halplotypecaller -vcf2- but I found there is a huge difference between both of them, for instance vcf1 was 4.4M and vcf2 was 140M.

When I tried to visualize both files I found that vcf2 has much more variants discovered than vc1.

Does anyone know why there is such difference? Does it have anything to do with gatk using known sites like dbsnp and freebayes don't?

I want to know because I was wondering if you can use FreeBayes instead of GATK in calling variants without affecting quality.

FreeBayes commandline:

freebayes -f human_g1k_v37.fasta out.bam > ~/out.vcf

GATK commandline:

```
java -jar gatk3.jar \
-T HaplotypeCaller \
-R human_g1k_v37.fasta \
-D dbsnp_vcf.vcf \
-o out.vcf \
-pairHMM VECTOR_LOGLESS_CACHING \
--emitRefConfidence GVCF \
--variant_index_type LINEAR \
--variant_index_parameter 128000 \
-A DepthPerAlleleBySample \
-stand_call_conf 30 \
-stand_emit_conf 10
```

Thanks in advance,
Shazly

Answer: I looks like you are assuming every line in vcf2 is a variant. However, you are using --emitRefConfidence GVCF which is asking HaplotypeCaller to make a call on every position, not just variants. If you are following Broad Best Practices, you need to also run GenotypeGVCFs in order to get the actual genotype calls. By the default GenotypeGVCFs will only output variant positions. This is the VCF you should compare to the Freebayes output.
Question: I have files of data that have Many of the columns but i need sequence of DNA ,I want to convert them to fasta then convert them to fastq or convert it directly to fastq , this files reach to gigabyte, what I must do ????
This is the data we have (tsv format) we need to convert it to FastQ [this is screenshot of it][1]

[1]: https://drive.google.com/file/d/0B_L8WR8veKL1WTV6Q1RQM2lpX0E/view

Answer: For those who may happen to reach this thread by way of search in future you can convert a fasta file to fastq format using `reformat.sh` from [BBMap][1] suite.

Please remember that the Q-scores created here are fake (example below sets Q-scores to 35 for all bases).

`reformat.sh in=test.fa out=fake.fq qfake=35`

[1]: https://sourceforge.net/projects/bbmap/
Question: I'm trying to give it a look at the phytozome database in biomart, but just as I want to explore a list of known loci, the program returns an error message related to biomaRt webservice. I have tried changing the Mart from phytozome to ensembl, and seems to be related to that. However, I don't see how that connects to the error given, do you have any hints to helping fix this?

This is my code:

library(biomaRt)
mart <- useMart(biomart = "phytozome_mart", dataset = "phytozome", host = "phytozome.jgi.doe.gov")
getBM(attributes = c("organism_name", "gene_name1"), filters = "gene_name_filter", values = "g400", mart = mart)
getBM(attributes = "organism_name", mart = mart)

Which returns error message:

>NULL
>Error in .processResults(postRes, mart = mart, sep = sep, fullXmlQuery = fullXmlQuery, :
The query to the BioMart webservice returned an invalid result.
The number of columns in the result table does not equal the number of attributes in the query.
Please report this on the support site at http://support.bioconductor.org

However, as mentoned above, I don't get this error if I try a similar approach with other databases, e.g. using Ensembl with pig data:

mart <- useMart(biomart = "ensembl", dataset = "sscrofa_gene_ensembl")
getBM(attributes = c("ensembl_gene_id"), mart = mart) %>% head

Should I really contact biomaRt support for this issue? Thank you

Answer: The solution here is to make sure you include `https://` in the host argument for Phytozome. If that isn't included the query is redirected back to the main Phytozome page rather than the BioMart service, R receives the home page HTML rather than a results table, and thus it fails.

library(biomaRt)
mart <- useMart(biomart = "phytozome_mart",
dataset = "phytozome",
host = "https://phytozome.jgi.doe.gov")

getBM(attributes = c("organism_name", "gene_name1"),
filters = "gene_name_filter",
values = "g400",
mart = mart)

#> [1] organism_name gene_name1
#> <0 rows> (or 0-length row.names)

The value "g400" doesn't seem like a valid gene name and we get an empty result, but if we use one that does exist we get:

getBM(attributes = c("organism_name", "gene_name1"),
filters = "gene_name_filter",
values = "82092",
mart = mart)

#> organism_name gene_name1
#> 1 Smoellendorffii 82092

I would not recommend running the version without anything supplied to the `values` and `filters` arguments. BioMart is not designed as a bulk data retrieval service and will probably fail (possibly without you realising by just omitting results). Better to register with Phytozome and download directly.

I'll update the biomaRt vignette to highlight that `https://` needs to be set for Phytozome.
Question: I followed the samtools/bcfutils/vcfutils pathway followed [here][1] to convert a set of human Hg19-aligned BAM files into a set of raw VCF files. I then got vcftools to filter down to just autosomal SNPs. These are really, really, really low-coverage genomes (they were enriched for NRY and/or mtDNA, and I am just trying to make use of the "leftovers")

Now I have the data I want, but I am trying to found out what of it is actually usable. I was wondering what are good filtering parameters for tossing/keeping human SNPs (or where can I find said parameters)? Thanks!

-Deven

[1]: http://ged.msu.edu/angus/tutorials-2012/snp_tutorial.html

Answer: This is what I use. I generally change them depending on the study. But more or less this is close to what everyone uses.

- MinDP (Minimum read depth): 5 (Indels) and 3 (SNPs)
- MaxDP (Maximum read depth): You have a low coverage data, so I would set it to 100. Normally it is 3 times the average coverage.
- BaseQualBias (Minimum p-value for baseQ bias): 0
- MinMQ (Minimum RMS mapping quality for SNPs): 20 or 30 (to be more stringent)
- Qual (Minimum value of QUAL field): 15 or 20

- StrandBias (Minimum p-value for strand bias): 0.0001
- EndDistBias (Minimum p-value for end distance bias): 0.0001
- MapQualBias (Minimum p-value for mapQ bias): 0
- VBD (Minimum Variant Distance Bias): 0 (More relevant to RNA-seq reads)

- GapWin (Window size for filtering adjacent gaps): 30 bp
- SnpGap (SNP within INT bp around a gap to be filtered): 20 bp

- SNPcluster (number of snps within a region): I usually drop all the snps if there are more than 3 snps within 10 bp.
Question: Hi all,

I am trying to index my bacterial genome using Hisat2 following the RNA-SEQ analysis workflow

I wonder which genome file format I should use for indexing.

I tried indexing both files;

1. cds_from_genomic.fna (named genome1)
2. _genomic.fna (named genome2)

*Using the command line ; hisat2-build genom1.fna indexed-genome*

as expected it generated 8 (indexed-genome.fna.ht2) files

indexed-genome.fna.1.ht2
indexed-genome.fna.2.ht2
indexed-genome.fna.3.ht2
indexed-genome.fna.4.ht2
indexed-genome.fna.5.ht2
indexed-genome.fna.6.ht2
indexed-genome.fna.7.ht2
indexed-genome.fna.8.ht2

Using the command line;

hisat2 -x indexed-genome -1 TF_paired.fq.gz -2 TR_paired.fq.gz -S T1-1_hits.sam

I got **T1-1_hits.sam** file, with size 0, along with this error ;

Could not locate a HISAT2 index corresponding to basename "indexed-genome"
Error: Encountered internal HISAT2 exception (#1)

So, did I choose the wrong genome-file.

Or I should use an alternative pipeline like using bowtie over Hisat2 for bacterial RNA-SEQ analysis?

Thank you in advance,

Answer: I think your basename might be "indexed-genome.fna".
Question: Hi,

I would like to ask if any of you know the callers of the pipelines that are stated on ICGC such as broad, embl and sanger? I checked the metadata however the website which directed me required id/password accession. These are pipelines that were used for calling mutations which produced VCF files.

Thank you

Best,
Tunc

Answer: Can you put the link you have been looking into? Usually the broad mutation caller is GATK toolkit and Mutect. Often VarScan2 is also used widely. But am assuming that this one you are looking into are the ICGC -TCGA project which is run by multiple sites in the world. I not sure if the publication is out yet. So [it shows only][1] the status of the calls made at all different contributing sites that are a part of this challenge so sequence WGS for large PAN-Cancer panel.

You can take a look at this [link][2] which has a somatic workflow for SNVs . This is of Sanger

I am not sure for the embl but Delly is one which is widely used but that is for SVs and not SNVs. I guess you have to do a bit more digging.

As I said you can find the below links for the different workflows that are being used by 3 major workflow managers , and I believe this is being circulated across different centers where samples are being analyzed.

- [The DKFZ/EMBL PanCancer Variant Calling Workflow][3]
- [DKFZ Dockered Workflows][4]
- [BroadWrapperWorkflow][5]

I think [this link][6] has most of the necessary information you are looking for.

[1]: http://pancancer.info/variant_calling.html
[2]: https://github.com/ICGC-TCGA-PanCancer/SeqWare-CGP-SomaticCore
[3]: https://github.com/ICGC-TCGA-PanCancer/DEWrapperWorkflow
[4]: https://github.com/ICGC-TCGA-PanCancer/dkfz_dockered_workflows
[5]: https://github.com/ICGC-TCGA-PanCancer/BroadWrapperWorkflow
[6]: https://github.com/ICGC-TCGA-PanCancer?page=3
Question: Of course, as a bioinformatician, I am aware of many large-scale open-source bioinformatics datasets, such as

- The ENCODE consortium (www.encodeproject.org, RNA-Seq, ChIP-Seq and so on),
- The Roadmap Epigenomics consortium (www.roadmapepigenomics.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The IHEC consortium (www.ihec-epigenomes.org, RNA-Seq, Chip-Seq, Bilsulfite-Seq),
- The TCGA/ICGC consortia (www.cancergenome.nih.gov, www.icgc.org, large-scale cancer data, DNA-Seq, RNA-Seq, etc.) and
- The LINCS consortium (www.lincscloud.org/l1000 , gene expression for more than a million of different perurbation experiments).

I am wondering, however, what other wonderful datasets, the are both large and open-source, are currently available. That might include things like RNA-Seq, Chip-Seq, Bisulfite-Seq, whole genome sequencing, WGAS, and many other assays (not necessarily NGS-related, though that is what I am mostly looking for).

Also things like the (neural) connectome of certain species (in any event large data) could be of interest.

There are quite some GEO datasets that at least partially fulfill these requirements, but most are simply having to few data samples in order to be interesting to me.

Your comments are greatly appreciated!

Answer: [BioGPS][1] is a good one for expression data.

[1]: http://biogps.org
Question: Hi All,
I have a question regarding trimming of adapters in NGS data. I have previously analyzed RNAseq data using brtim (https://www.sciencedirect.com/science/article/pii/S0888754311001339) without using adapter sequence. I am now analyzing smallRNAseq data, and I am using adapter sequence as adapter.fa using BBduk.sh from BBmap tools. Could you please clarify under what circumstances one would need to provide/know adapter sequence and when it is not necessary?

Answer: This is a good question for the beginners who are starting with the NGS analysis.

> In common short read sequencing, the DNA insert (original molecule to
> be sequenced) is downstream from the read primer, meaning that the 5'
> adapters will not appear in the sequenced read. But, if the fragment
> is shorter than the number of bases sequenced, one will sequence into
> the 3' adapter. To make it clear: In Illumina sequencing, adapter
> sequences will only occur at the 3' end of the read and only if the
> DNA insert is shorter than the number of sequencing cycles (see
> picture below)!

This is true in case of smallrna sequencing like microRNAs and piRNAs which are between 16-33 nt long. That means even with 50bp single end sequencing, 3' adapter sequences will occur. This is generally not the case with RNAseq where the gene/transcript fragments are way longer than the sequencing length and adapter sequences do not accur.

![enter image description here][1]

For details: https://www.ecseq.com/support/ngs/trimming-adapter-sequences-is-it-necessary

Also, if you are not sure which adapters are used for the sequencing, you can use an adapter predicting software like DNApi that will predict the most common adapter sequence which can be used by the tools like cutadapt or STAR before the mapping.
![enter image description here][2]
Paper: Tsuji, J., & Weng, Z. (2016). [DNApi: a de novo adapter prediction algorithm for small RNA sequencing data][3]. PloS one, 11(10), e0164228.

link: https://github.com/jnktsj/DNApi

I hope this helps.

[1]: https://www.ecseq.com/support/ngs/img/fragmentsize.png
[2]: http://journals.plos.org/plosone/article/figure/image?size=large&id=10.1371/journal.pone.0164228.g001
[3]: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0164228
Question: I'm searching for the most comprehensive source of chromatin interactions to support enhancer targets (data such like Hi-C, ChIA-PET, IM-PET, 5C, 4C, 3C, etc.) and my question according that could be split into two following parts:

1. On page [http://epigenie.com/epigenetic-tools-and-databases/][1] I've found out links to two databases, but the reference to 4D Genome doesn't work. **Does anybody know whether this database
[http://4dgenome.int-med.uiowa.edu][2] was moved somewhere else or maybe its only internally available for U Iowa community?** (I'm asking about probably exclusive access, because on Tan's lab page [http://www.healthcare.uiowa.edu/labs/tan/software.html][3] only the 4DGenome has the trailing part of the link *int-med.uiowa.edu* (presumably this **int-med** strands for internal/general medicine or maybe in the sense of internal/intranet access), whereas all other services like IM-PET has its link starting with *http://www.healthcare.uiowa.edu/labs/tan/*.

2. **Maybe you know even better database for chromatin interactions covering as much as possible cell lines/tissue types for GRCh38 assembly or any other whose coordinates could be easily converted, e.g. by *liftOver*?** I need data only for human and mouse, but I would be grateful also for solely human data.

Thank you in advance for your answer.

[1]: http://epigenie.com/epigenetic-tools-and-databases/
[2]: http://4dgenome.int-med.uiowa.edu
[3]: http://www.healthcare.uiowa.edu/labs/tan/software.html

Answer: This problem has been finally solved with help from prof. Kai Tan.
It occurred that his lab has been moved from University of Iowa to Children's Hopspital of Philadelphia, so has been this service.
Now, the working link is:
**[http://4dgenome.research.chop.edu][1]**

[1]: http://4dgenome.research.chop.edu
Question: I've produced a set of about 400 of GVCF files with gatk HaplotypeCaller, with the `-ERC GVCF` option. I'd now like to combine them for downstream genotyping and variant recalibration. I believe I can combine with gatk CombineGVCFs.

gatk CombineGVCFs \
-R reference.fasta \
--variant sample1.g.vcf.gz \
--variant sample2.g.vcf.gz \
-O cohort.g.vcf.gz

But what I don't know, is how to input all my 400 GVCF files into CombineGVCFs. I've heard this can be done with the `--arguments_file` option, but I don't know how to build such a file?

Any help gratefully received!

Answer: create a file with the `.list` suffix containing the path to your vcf

find /path/to/dir -type f -name "*.vcf.gz" > input.list

and then use this file with the '-V' argument

(...) --variant input.list
Question: can anyone tell me the difference between lncrna as targets of mirna and lncrna as target mimics of mirna?

Answer: In the former case, the miRNA is being targeted toward the lncRNA to induce RNAi of it. In the latter case, the lncRNA is acting more like a sponge, but complementing the miRNA and thereby preventing it from binding to its normal target. In reality, both could and likely do occur at the same time. It would be unsurprising to find that a lncRNA has some notable function in and of itself (thus making it an interesting miRNA target), but that the same transcript could also act as a sponge for miRNAs, preventing them from interacting with other targets. Biology likes to be messy like that.
Question: I'm pretty noob at coding and I finally have a problem where I absolutely want to automate.

I'm interested in a 50-or-so amino acid motif that occurs multiple times in the same protein. I can do a multiple sequence alignment of the 4 times it occurs in my favorite protein that reveals sequence conservation. It has the form:

(6-residue conserved motif)-(variable loop of ~6-20 AAs)-(12-residue conserved motif)-(short variable 5-7 AA loop)-(10-residue conserved motif)

Based on the conserved features I can see in the alignment, I can identify this motif in proteins in the NCBI database. My main goal is to be able to supply an input sequence and have the script identify all instances of the motif. My first question is: do I need to build my own hidden Markov model or are there likely tools out there that would do the trick?

I can do a multiple sequence alignment with my 4 known motifs and an unknown sequence that is expected to have the motif. This seems to help the alignment software find the pattern in the unknown sequence, if it exists. However, it only finds the best match in the unknown sequence, and I have to remove that motif from the input to then find the 2nd best-matching motif, and so on. This obviously takes way too long, but I realized I definitely want to somehow incorporate the information from a multiple alignment of motifs into the way I search for new ones. This is especially true because using my initial alignment to search for the motif often finds ones that are notably different from the ones in my favorite protein, but still have the same overall patter (e.g. hydrophobic-serine-hydrophobic-proline becomes hydrophobic-threonine-hydrophobic-proline where the hydrophobics are different and serine becomes threonine, but the pattern is the same and in in the same position in the 50AA motif). So upon finding new ones, I can align them to the ones from my favorite protein, resulting in a new alignment that better represents the possible sequence space of this motif. In my mind, as I iteratively add to my alignment, it will become better and better at identifying different versions of this motif. What would be the best way to continually update the definition of the motif that I'm searching for using the ones I've identified?

Answer: Some papers to consider:

- https://academic.oup.com/bioinformatics/article/24/6/807/194276
- https://academic.oup.com/bib/article/15/4/582/408173
- https://academic.oup.com/nar/article/33/suppl_2/W239/2505540

And websites:

- https://toolkit.tuebingen.mpg.de/tools/hhrepid
- https://toolkit.tuebingen.mpg.de/tools/repper

A tool like [**HMMer**][1] can be used to create a custom hidden Markov model, which will find multiple repeat occurrences in a single pass. You may have to run it more than once to identify all repeats, but it is reasonably fast.

[1]: http://hmmer.org/
Question: I'am trying to open index.html from JBrowse after running ./setup.sh.
I have just loaded sample_data of volvox.
But some of the tracks are not loading file.
Showing "Error : HTTP 200 when fetching file"

Please help me in this
Thank You.!

Screenshot : https://photos.app.goo.gl/NQsa5SJC5w7APfHk8

Answer: The setup data should work, however, it looks like you are trying to run it on your local client machine using the file:/// protocol in the browser?

If that is true, try using the JBrowse-1.16.1-desktop-linux-x64.zip from the following site, it is likely more appropriate to your needs.

http://jbrowse.org/blog/

I use JBrowse on a server only, but the local version is attractive for some colleagues.
Question: I would need to retrieve the normalisation (and maybe the batch correction method) used to produced the pan-Cancer Atlas mRNA expression matrix (file called 'EBPlusPlusAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.tsv' found [here][1]).

Starting from the raw read counts obtained from the GDC and the same gene panel, I tried FPKM and FPKM-UQ normalisation as described [here][2], but the expression values obtained do not fall at all in the same range as in the pan-Cancer mRNA matrix. Maybe that would suggest a cross-sample batch correction.

My goal is, starting from raw read counts, to normalise expression data from new samples together with the pan-Cancer mRNA data, in order to get a unified expression matrix and to be able to compare apples to apples basically.

Any information or alternative method would be greatly appreciated.

[1]: https://gdc.cancer.gov/about-data/publications/PanCanStemness-2018
[2]: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/#mrna-expression-workflow

Answer: My guess (and it is only a guess), given the name of the file, that this is build from the RSEM quantification results that are present in the Broad Institute's [Firehose porta][1]l, rather than from read counts.

RSEM use an EM algorithm to build isofrom expression values. A length-weighted sum of these values is then used to create gene expression values.

The firehose documentation states that these are normalised like so:

> RSEM expression estimates are normalized to set the upper quartile
> count at 1000 for gene level and 300 for isoform level estimates.

Please note that this is definately NOT a batch correction and that [batch effects have been shown to be a serious problem with PanCancer analyses][2] (although this is at the level of somatic variants)

[1]: https://gdac.broadinstitute.org/
[2]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5387285/