cannin/biostars_qa huggingface dataset:
Question: There are number of Bioconductor packages published. If anyone have some problem, how can they really contact the authors of the package, to get the help or to ask about their doubts?

After getting stuck a number of times, I google the problem, and if I don't find answer, I write email to the authors. Among the emails I wrote in past, I received only one response for my problem, that I was facing with [qusage][1] package. I am not mentioning about the package for which I didn't get any response from authors.

I know that one way could be to to post query on [Bioconductor Support][2]. Can there be a online forum or service, by which, the authors of the package also get an online page for that package, to ask question, may be at the end of the page of their package on Bioconductor website. I believe that author of the package would be able to able to answer more clearly to the questions of users.

Just creating a package, publishing it and forgetting is not good!

I have also posted the same on [Bioconductor Support][3] website.

[1]: http://bioconductor.org/packages/release/bioc/html/qusage.html
[2]: https://support.bioconductor.org/
[3]: https://support.bioconductor.org/p/68914/

Answer: As Dan mentioned in [his answer][1] on [the bioconductor support site][2], just post the question on the support site and tag it with the name of the package. They have implemented [maintainer integration with the support site][3] which tries very hard to tickle the listed maintainer of the package via email whenever a question is posted that is tagged with a package they maintain.

If the maintainer doesn't tune in to the website to provide an answer, then you're already talking to an audience that would be your 2nd best bet to get an answer from.

[1]: https://support.bioconductor.org/p/68914/#68915
[2]: https://support.bioconductor.org
[3]: https://stat.ethz.ch/pipermail/bioc-devel/2015-May/007532.html
Question: I have a fastq file and a fasta file with one sequence. I want to plot a figure to show all the reads in fastq file aligned to the fasta file as shown in the figure below. I want to plot for both plus and minus strands. How can I do this (preferably using R)? Thank you for your help.
Plot I want:
https://ibb.co/eV7EDx

Answer: Just a suggestion:

1. Align your fastq files against the fasta file, perhaps using bwa or bowtie2 outside of R.
2. Read the resulting BAM files into R using the GenomicAlignments Bioconductor package.
3. Visualize data using the Gviz Bioconductor package.
Question: What are chimeric reads? Is there any effective method to distinguish chimeric reads from normal ones? Can someone point me to a survey or reviews on that topic? Thank you a lot!

Zhen

Answer: Chimeric reads occur when one sequencing read aligns to two distinct portions of the genome with little or no overlap. Chimeric reads are indicative of structural variation. Chimeric reads are also called split reads.

After aligning with bwa mem, chimeric reads will have an SA tag as described on page 7 of the [SAM format specification][1]. To find them all you have to do is extract them using grep.

samtools view my_alignment.bam | grep 'SA:' | less

For reference:

- [Identification of genomic indels and structural variations using split reads][2]
- [Genome structural variation discovery and genotyping][3]

[1]: http://samtools.github.io/hts-specs/SAMv1.pdf
[2]: http://www.biomedcentral.com/1471-2164/12/375/
[3]: http://www.nature.com/nrg/journal/v12/n5/abs/nrg2958.html
Question: How can I filter delly output to only deletions and ones that PASSED the quality check? I want to use the delly output from delly call.

Answer: [Bcftools][1] works best for filtering VCF/BCF files. For instance

`bcftools query -f '%CHROM\t%POS\t%INFO/END\t%ID\t%INFO/SVTYPE\t%QUAL\t%FILTER\n' delly.bcf`

generates a table with all the information you want. Alternatively, you can also filter the VCF/BCF directly using:

`bcftools view -i 'FILTER="PASS" & INFO/SVTYPE="DEL"' delly.bcf`

[1]: http://samtools.github.io/bcftools/
Question: I want to visualize the variants I have in my VCF (or any tab delimited) file.

Specifically, I want to look at them in relation to annotation tracks such as Cons Score, Encode TFBS binding site, exon or no (and if so what amino acid) and many more.

This can be done through UCSC: https://genome.ucsc.edu/cgi-bin/hgTracks?hgsid=389193029_56wIWP5VNLUDfDyneR3PpEQOJdyN&hgt.psOutput=on

possibly through IGV? And probably through many other tools.

*I am writing to ask the community which tools make the best-quality and most attractive graphs for visualizing NGS data?*

I would be looking at something between a few hundred bases up to a whole gene (100s to 100000s of bases) so I do not need something like Circos.

Suggestions? Thank you.

Answer: In addition to Gviz, as mentioned by https://www.biostars.org/u/6141/, there's also ggbio. I find both packages pretty useful.
Question: I am aligning RNA-Seq data to human genome reference hg38. How can I build the index using hg38.fa

Thanks in advance

Answer: Bowtie
--------

Building a new index

Bowtie2
-------

The bowtie2-build indexer
Question: Hi All,
I need some help interpreting my DEseq2 results. Here is a bit of my DEseq2 code:

dds$cond <- relevel(dds$cond, ref = "C")
dds = DESeq(dds, fitType = "parametric")
my.names<-resultsNames(dds)
my.names
[1] "Intercept" "cond_T_vs_C"

Then I did:

res<- results(dds, contrast=list( "cond_T_vs_C"), test="Wald")

Now looking at the result, for this particular gene I have:

baseMean log2FoldChange lfcSE stat pvalue padj
Seq_483975_2034 1995.245 -1.354257 0.4045764 -3.347346 0.000815893 0.005723768

Can someone please explain me whether this gene is down regulated in Control (C) or in Treatment (T)? Thanks

Answer: Use the `plotCounts()` function to inspect the gene of interest, you'll be able to see then if the gene is up- or down-regulated.
Question: Hello,

If there is a SNP that surpassed GWAS significance, but it turned out to be triallelic, should we need to exclude this variant? I merged the two datasets, cases and controls, but no issues in the process to analyze in plink. But dbSNP says it is triallelic. In the map files for the two datasets, the SNP's alleles are both A/C. But dbSNP says it is A/C/T (T is ancestral).

Thank you in advance!

Answer: Do you know how the genotypes were initially measured? If the third allele could be expected to be systematically mistaken for one of the first two, you should adjust your interpretation of the result accordingly. But if the third allele legitimately doesn’t seem to be in your dataset and is generally rare in the population you’re studying, I wouldn’t worry about it.
Question: Hi,

I've been trying to run ABySS on PE data (~90GB) from a plant genomic sample with a k-mer size of 70 (k=70) on an HPC system. For memory concerns, I would like to use paired de Bruijn approach; however, I am not sure how to pick appropriate individual k-mer size (K) for a k-mer pair span (k) of 70 (or, any K&k values to get a close result to single k=70). I suppose larger gap between pairs means less memory, but is there a way to estimate K for a reasonable gap size?

[For those who may have some more spare time;

In fact, I don't think I get the idea with the k-mer pairs at all. How can the gap be larger than the K value (in the readme, for E.coli there is an example with K=16 and k=64, making the size of the gap 32). In this mode, does the graph again looks for K-1 overlaps? Any explanations or sources, from where I can read more, will be of great help.]

Many thanks in advance, any help is greatly appreciated.

Best,
Ani

Answer: Hi Ani,

I don't know any easy way to choose optimal values for K and k; probably the most practical thing is to do is assemble with a range of values and see what gives the best results.

Here is an explanation of the paired de Bruijn graph idea, though, which may help with your intuition. Consider an example k-mer size of 8, e.g:

ACGTACGT

With ABySS paired de Bruijn graph, the parameters K=3 (individual k-mer size) and k=8 (k-mer pair span) would give you a paired k-mer like this:

ACGNNCGT

where the N's are "wildcard" positions (in other words they can match any base). The intuition of the paired de Bruijn is that, in most cases, the standard 8-mer and the paired-kmer will match the same places in the genome, so it doesn't matter that we are "throwing away" those bases in the middle.

By that reasoning, 'k' (k-mer pair span) is generally a more important parameter than 'K' (individual k-mer size). So when you are doing your parameter sweep, it it probably best to fix 'K' at some value (e.g. 32) and assemble with a range of 'k' values. Then, when you find a good 'k', you can try adjusting 'K'.

The paired de Bruijn graph is contructed by creating a node for each paired k-mer and an edge between that nodes that have *paired* overlaps.

The paired de Bruijn graph idea was described in:

Medvedev, Paul, et al. "Paired de bruijn graphs: a novel approach for incorporating mate pair information into genome assemblers." *Journal of Computational Biology* 18.11 (2011): 1625-1634. URL: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3216098/

Also, this book has some nice diagrams explaining the idea:

Jones, Neil C., and Pavel Pevzner. *An introduction to bioinformatics algorithms*. MIT press, 2004. Website: http://bioinformaticsalgorithms.com/index.htm

Pavel Pevzner et al. also had an online course that explained the paired de Bruijn, but I can't seem to find it anymore. Maybe you will have better luck.
Question: When selecting a commercial primer for a particular organism, is it necessary to validate that the primers can be amplified amongst the genomic DNA of the host organism? Wouldn’t commercial primers typically be validated in this way already? Also, even if I were to design my primers in silico, would I need to validate amplification amongst genomic DNA?

Answer: Yes it is good practice to validate primers on one or a few samples before undertaking your experiment. Things to check are whether the PCR product is the expected size. While the primers are likely tested by the supplier, you want to be certain that it is compatible with the elution buffer all of the other PCR reagents. If you are designing your own primers, you need to check that only 1 product is generated. What is your application? Is it RT-qPCR or something else?
Question: Hi all! I wasnt successful trying to find files with annotated transposable elements on chromosomes from Arabidopsis. I am also interested in their specific locations on chromosomes. Any help will be appreciated :)

Answer: Your best bet for Arabidopsis genome annotation is TAIR, and the file you are looking for is [this one (FTP link)](ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_transposable_elements/TAIR10_Transposable_Elements.txt).

These data are based on the work of Buisine, Quesneville and Colot (Genomics 2008), for more information see: http://urgi.versailles.inra.fr/Data/Transposable-elements/Arabidopsis
Question: Hello everyone,

I have a set of differentially expressed genes in an "A" tissue when compared with other tissues, most of these genes being typical protein coding genes, and the rest of them being transcription factors. I have obtained this set from RNA-seq data from several samples of this "A" tissue and the rest of the tissues. This far, I can assume these typical protein coding genes are regulated by the transcription factors.

However, my problem is, is there any way or tool to know which genes are regulated by each transcription factor without using any database or experiment, but only my data?

Thank you very much for your help and kind regards.

Answer: If we rule out using any database (I don't know why we'd do that, using a database of TF binding sites would be convenient), then one possibility would be to look at covariation between the various DE genes (e.g., with WGCNA). The assumption would be that a TF driving a bunch of genes to be DE should be near the center of a module containing them and possibly be its eigengene.
Question: Hi all,

I'm currently struggling with structural variant calling.

I would like to call those structural variants for all samples in my population of haploid fungi (around 100).

My aim is to cross the results of all samples to get absence/presence polymorphisms of all called variants.

Would anyone know a tool that is able to call variants on all samples, and then cross the results to have only one file with SVs presence/absence for all samples?

Thanks in advance for your help

Guillaume

EDIT : I've eventually been using Delly2, which can handle multi-sample calling

Answer: If you are planning to look for Indels mainly, then [Freebayes][1] performs well for population data. But if you also interested in other SVs too, then you would need one of the [SV-caller pipelines][2] seperately on each of the populations and then look for those which are common in all, since calling it in one go for all populations would just go to majority voting through read-support - not sure if any of these are population based though

[1]: https://github.com/ekg/freebayes
[2]: https://omictools.com/structural-variant-detection-category
Question: I am absolutely new to Structural Bioinformatics (only started last week). I am working on the secondary structure assignment/prediction (actually I am not sure) of proteins using ML.

When it comes to protein analysis, most of the books or research papers talk about protein structure **Prediction**. Only a handful of books or articles talk about protein structure **Assignments**.

Can anyone tell me what is the basic difference between Assignment and Prediction?

Answer: Assignments are made from known 3D structures, and are fairly straightforward. Predictions are made without 3D structures. They are usually based on residue conservation profiles that are extracted from multiple sequence alignments.
Question: Hello,

Does anyone know where I can find a list of Ensembl transcript IDs (e.g., ENST0000011111) associated with protein-coding genes?

I've run my mRNA-Seq data through kallisto and am plotting PCA on the pseudocounts, which are associated with Ensembl Transcript IDs. A labmate wonders if my PCA would cluster better if I only used the rows of counts corresponding to Ensembl Transcript IDs associated with mRNA (since, presumably, all of the other rows would be 0, since the sequencing was only done on mRNA, not ncRNA).

I'm not sure to what extent removing rows of 0s would change my PCA, but I thought I'd give it a try.

Thanks for your help!

Answer: Go to:

http://www.ensembl.org/biomart/martview/21893b6a5bf91112c9a786f296b9200d

Choose database > Ensembl Gene 83 > Homo sapiens

Attributes > Gene > click on Ensembl Protein ID (or any other protein IDs such as refseq, etc.)

go to results and download
Question: Dear all,

Do you have any idea how to easy extract contigs from fasta file wich contains specific sequence?

For example:

My sequence:

ACCGTACCC

My FASTA:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
>c1044
ATGAATAAAATAATTTTGTATCATAAATCGAGCTGTTAATTATT
>c1044
TTCATATTTGTAGCTAAGCAGAGGCGAAGCGTTCTTGTATCG
```

My output:

```
>c1042
ACCGTACCC
>c1043
GCTACAGTTGAAAGGGGACCGTACCC
```

Thank you so much for any ideas and help.

Answer: https://gist.github.com/mdshw5/c6b9f7c08d36b3a3c5b0

In this case I've changed the duplicate defline since pyfaidx requires unique sequence ids. You could mangle all the key names by passing your own `key_function` if you like.
Question: I invested some time into learning R and would prefer not to learn Python or Ruby just for convenience of parsing XML.

I am not a pro R programmer so I need a higher level library then XML package of R.
getNodeSet and XPath is very hard to use.

Answer: As for biophyton and bioperl counterpart in R you can use Bioconductor, however it may not have all the functionality you find in those libraries, and bioconductor is the project, you would need to search in biconductor's website for specific packages.
Question: I am using plink to ld-pruning SNPs. But I wonder why plink always produce the same pruning results, even when I add a random seed.

Let's suppose, SNPs A and B are in LD. They can be also in LD with other SNPs. The SNPs in LD with A and the SNPs with B might be different. Thus, the first choice of A or B might produce different pruning SNP set. So, LD-pruning should produce a random result, really?

Here, I worry that Plink always produces a same result, it will produce a bias? Thank you for comments.

Answer: PLINK keeps the SNP with higher minor allele frequency, if the two MAFs are not identical. When there is a tie, the first SNP is kept.
Question: Background: The CFTR gene has a Poly TG/Poly T region in intron 8 that has [varying clinical][1] consequences.
I'm using targeted NGS for this region (150bp, paried end) and I'm trying to call the PolyTG/PolyT genotypes from fastq files. My first attempt was bbmap run with default parameters and a reference including (TG)9-12(T)5,7,9:

bbmap.sh ref=$ref in1=$R1 in2=$R2 scafstats=$out

39/40 samples could be correctly called based on the bbmap stats (header shown below). Heterozygous calls were made from ref_name's with the highest unambiguous reads. Homozygotes could be figured based on the ratio of the top two.

#sample	ref_name	per_unambiguousReads	unambiguousMB	per_ambiguousReads	ambiguousMB	unambiguousReads	ambiguousReads	assignedReads	assignedBases

I expanded the reference to catch rare genotypes (e.g. T6 and TG13). The expanded reference has (TG)8-13(T)2-10. Running bbmap with this reference results in many incorrect calls.

I'm wondering if there are parameters I should be tweaking in bbmap or if there is a better solution to calling the Poly TG/Poly T genotypes that I'm not thinking of.

Thanks,
Bob

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2699602/

Answer: Just to wrap this up. I found this [paper][1]. I contacted Dr. Pagin and he explained their script used grep on fastqs and counted different alleles to make a call per sample. I used that idea in this [script][2]. If anyone uses this script in the future, the HET_CUTOFF (line 554) will likely vary depending on the panel you use to sequence CFTR; adjust according to your data.

[1]: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0149426
[2]: https://github.com/rjsicko/scripts/blob/master/polyTG_table_v1.0.sh
Question: Samtools can be used to select reads above certain mapping quality.

samtools view -h -b -q 30 aligned.bam -o above.mapQ30.bam

**But, how to select a read below certain mapping quality - all aligned reads below mapQ 30?**

I know it can be done using awk. But, the pipeline gets lengthy and time consuming when first need to convert bam to sam - separate header - use awk for mapQ below 30 - add header - sam file - convert to bam.
Really, its taking a lots of time.

Thanks,

Answer: Use SAMtools with the -U flag to output reads that don't match your filter, e.g.

> samtools view -h -b -q 30 -U below_q30.bam aligned.bam
Question: I need to upload several terabytes of files to SRA. I read about something called "ascp" that works with the Aspera server NCBI uses but I can't seem to download this software without a user/pass.

Does anyone have this software?

Answer: I had this same issue (and asked myself why anyone would use this over rsync), and then found the answer. It is extremely non-intuitive, but you must download the Aspera Connect browser plugin, and the CLI utility `ascp` is included in that package. If you still have issues, I have it in my Dropbox.

**Update**: Aspera now maintains this shell script as a separate file. You should now be able to do

sh <(curl -s http://demo.asperasoft.com/ascp-install-3.5.4.102989-linux-64.sh)
Question: Dear all,
I am analyzing my single cell data using Scanpy package and as you may knwo it stores my data in anndata object which is a panda like data frame. I am interested to know which of my cells express the gene x. Does anyone has an idea how would I be able to access that information? I am not really familiar with data structures in Python and specially panda

Answer: Ok I found my answer :D It was not that complicated

adata[:, adata.var_names.str.match('genenames')]
Question: I want to annotate gene names for a segmented file by finding overlap with it with a reference file. Would like to know how can it be done using GenomicRanges Package or any other bioconductor packages.

`Segmented_file`

```
Chromosome Start End
1 5930000 11730000
1 16850000 18010000
```

`reference_file`

```
Chr Start End Gene
1 5930500 6230500 SPSB1
1 6930500 7340500 SPSB2
1 16854500 16950000 TAS1R1
1 17810032 17910064 ENO1
```

Expected results

```
Chromosome Start End Gene
1 5930000 11730000 SPSB1,SPSB2
1 16850000 18010000 TAS1R1,ENO1
```

Answer: Yes, you can use the GenomicRanges library in R to do this. Here is an example:

```r
library(GenomicRanges)
# generate 10 random segments
regions.gr <- GRanges(seqnames=Rle(rep(1,10)), ranges=IRanges(start=sample(1:1000,10), width=100), strand=Rle(rep("*",10)))

# generate 5 random genes
fiveGeneNames <- paste(sample(letters,5), sample(1:100,5), sep="")
genes.gr <- GRanges(seqnames=Rle(rep(1,5)), ranges=IRanges(start=sample(1:1000,5), width=100, names=fiveGeneNames), strand=Rle(sample(c("+","-"),5, replace=T)))

# what regions overlap what genes?
overlapGenes <- findOverlaps(regions.gr, genes.gr)

# Return any genes with an overlap.
# Convert the resulting "Hits" object to a data frame
# and use it as an index
overlapGenes.df <- as.data.frame(overlapGenes)
names(genes.gr[overlapGenes.df$subjectHits])

# extract the regions that hit genes
regionsWithHits <- as.data.frame(regions.gr[overlapGenes.df$queryHits])
# add the names of the genes
regionsWithHits$genes <- names(genes.gr)[overlapGenes.df$subjectHits]

# Some segments might not overlap any gene.
# Find the distance to the nearest gene
distanceToNearest(regions.gr, genes.gr)

# if your genes and segments were in a bed file you could easily import them
library(rtracklayer)
regions <- import.bed("myRegions.bed",asRangedData=FALSE)
genes <- import.bed("myGenes.bed",asRangedData=FALSE)
```
Question: I can make box plot using the base R plot ,but I want to do the same with ggplot2 , its not as simple as making it in normal boxplot function ,

This is my data frame

Gene `7256_Mono` `7653_Mono` `6792_Mono` `6792_HSC` `7653_HSC` `7256_HSC`
DNMT3A 9.8143600 5.5681400 3.32429000 7.7649600 7.8304600 21.9411000
DNMT3B 0.1053070 0.0826214 0.20910300 2.3988700 3.3696100 3.9184700

So this is my subset of data set how do I coerce this dataframe into ggplot data format , i have read tutorials but it quite confusing ,I have to make factors which is used in the ggplot , so how do i make factors out of this data sets so that I can use it for ggplot and make box plot out of it. Any help and suggestion would be highly appreciated .

Answer: Hi krushnach80, not sure how you want to group the variables, by but below is my solution.

> head(dat)
Gene X.7256_Mono. X.7653_Mono. X.6792_Mono. X.6792_HSC. X.7653_HSC. X.7256_HSC.
1 DNMT3A 9.814360 5.5681400 3.324290 7.76496 7.83046 21.94110
2 DNMT3B 0.105307 0.0826214 0.209103 2.39887 3.36961 3.91847

## Libraries... I think there is one library for all these now.
library(tidyr)
library(dplyr)
library(ggplot2)

## Format the data as 'long' using dplyr %>% and tidyr gather, and omitting Gene.
mdat = dat %>% gather(key = sample, value = val, -Gene)
head(mdat)
Gene sample val
1 DNMT3A X.7256_Mono. 9.8143600
2 DNMT3B X.7256_Mono. 0.1053070
3 DNMT3A X.7653_Mono. 5.5681400
4 DNMT3B X.7653_Mono. 0.0826214
5 DNMT3A X.6792_Mono. 3.3242900
6 DNMT3B X.6792_Mono. 0.2091030

## Make boxplot with scatter overlay because boxplots can be misleading. Grouping factor is gene. Also try beeswarm-type plots.
ggplot(mdat, aes(Gene, val)) + geom_boxplot(colour = "grey", fill = "light grey") + theme_bw() + geom_jitter(width = 0.2, alpha = .9, shape = 2)