cannin/biostars_qa huggingface dataset:
Question: Dear everyone,

I have been asked to comment on an experiment design that involves WES or WGS of cancer cell lines which lack matched normal.

I know that this design is far from ideal, but I was wondering if there are people who have already stream-lined it (to the extent possible of course). I can think of a couple of common variation filters to keep out the germline variants, but I am not sure how to go beyond.

I am sorry if this is a duplicate question; I couldn't find a related hit.

Thank you,

Noushin

Answer: Since you are asking about cell lines, maybe these have already been sequenced. There are bigger studies I'm aware of:

1. Cancer Cell Line Encyclopedia (CCLE), see http://www.ncbi.nlm.nih.gov/pubmed/22460905
Browse and download the data: http://www.broadinstitute.org/ccle/home
2. NCI-60 cell line, see e.g. here: http://www.cbioportal.org/public-portal/study.do?cancer_study_id=cellline_nci60

Then I would check against COSMIC for known somatic mutations. And maybe the cell line in question even has data in COSMIC, see By Sample at http://cancer.sanger.ac.uk/cancergenome/projects/cosmic/
Question: Hi all,

I am trying to subset my extremely large vcf file of around 11,000 samples based on the string "pathogenic". I tried the following:

bcftools view -i 'CSQ="Pathogenic"' -o PathOnlyClinVar.vcf Merged.vcf.gz

which I believe should have worked, but it only returned the header without any other information. I mainly want to safely subset my vcf file with ClinVar information that contains the string "Pathogenic".

Thanks!

Answer: https://samtools.github.io/bcftools/bcftools.html#expressions

> regex operators "\~" and its negation "!~". The expressions are case sensitive unless "/i" is added.

bcftools view -i 'CSQ ~ "Pathogenic"' -o PathOnlyClinVar.vcf Merged.vcf.gz
Question: Hi,

I have analysis with Qiime of 16s reads, and I want to make a pie chart with Krona but do not know which file from Qiime to use for Krona?

Answer: > Hi, hope this answer is not too late. :)

## convert biom to tsv format

biom convert -i single_sample.biom -o single_sample.tsv --to-tsv --table-type "OTU table" --header-key taxonomy

## remove the first two rows with comment

# take a look at single_sample.tsv
$ head single_sample.tsv
# Constructed from biom file
#OTU ID	H2O	taxonomy
346085	140.0	Bacteria; Proteobacteria; Alphaproteobacteria; Caulobacterales; Caulobacteraceae; Brevundimonas; Brevundimonas_bullata
10298	2.0	Bacteria; Proteobacteria; Gammaproteobacteria; Enterobacteriales; Enterobacteriaceae; Photorhabdus; Photorhabdus_temperata
122823	3.0	Bacteria; Proteobacteria; Betaproteobacteria; Burkholderiales; Oxalobacteraceae; Massilia; EF516371_s
130468	2.0	Bacteria; Proteobacteria; Alphaproteobacteria; Sphingomonadales; Sphingomonadaceae; Sphingopyxis; Sphingopyxis_witflariensis
139977	38.0	Bacteria; Proteobacteria; Alphaproteobacteria; Sphingomonadales; Erythrobacteraceae; Erythrobacter; Erythrobacter_flavus
121751	2.0	Bacteria; Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Catonella; JX096343_s
96934	10.0	Bacteria; Proteobacteria; Gammaproteobacteria; Pseudomonadales; Pseudomonadaceae; Pseudomonas; Pseudomonas_guguanensis
95181	4.0	Bacteria; Proteobacteria; Gammaproteobacteria; Pseudomonadales; Moraxellaceae; Acinetobacter; Acinetobacter_radioresistens

# remove the comments
$ egrep -v "^#" single_sample.tsv > sample_no_comment.tsv

## remove the first column

$ cut sample_no_comment.tsv -f 2- > sample.tsv
$ head -n2 sample.tsv
140.0	Bacteria; Proteobacteria; Alphaproteobacteria; Caulobacterales; Caulobacteraceae; Brevundimonas; Brevundimonas_bullata
2.0	Bacteria; Proteobacteria; Gammaproteobacteria; Enterobacteriales; Enterobacteriaceae; Photorhabdus; Photorhabdus_temperata

## replace the "; " with "\t"

$ sed -i 's/;\s*/\t/g' sample.tsv
$ head -n2 sample.tsv
140.0	Bacteria	Proteobacteria	Alphaproteobacteria	Caulobacterales	Caulobacteraceae	Brevundimonas	Brevundimonas_bullata
2.0	Bacteria	Proteobacteria	Gammaproteobacteria	Enterobacteriales Enterobacteriaceae	Photorhabdus	Photorhabdus_temperata

## now we can make a pie chart with ktImportText(One of the Krona tools)

$ ktImportText -n enjoy_your_life -o sample.krona.html sample.tsv

## open sample.krona.html with your web browser, you'll see this
![Krona_snapshot_enjoy_your_life][1]

## all in one command(optional)

$ awk -F '\t|;' 'BEGIN{OFS="\t"} FNR > 2 {$1=""; print $0}' single_sample.tsv | cut -f 2- > sample.tsv

[1]: https://preview.ibb.co/hLhkoQ/Krona_snapshot_enjoy_your_life.png
Question: Hi,

I was wondering if there is a possibility of marking genes of interest on the heatmap via `ComplexHeatmap` R package. Instead of specifying row numbers in` anno_mark` function, can we specify the gene list and mark on the heatmap? For long gene list, it's bit challenging to identify the rownames of each gene and specify them. Thank you.

## Genes of interest
genes_to_show = c('Gene_A', 'Gene_F', 'Gene_L', 'Gene_R', 'Gene_GG', 'Gene_ZZ', 'Gene_Test')

## Mark rownames
ha = rowAnnotation(foo = anno_mark(at = c(1:4, 11, 28, 44,:47),
labels = rownames(m)) ## Mention gene list here "genes_to_show" instead of specifying row numbers

## This would work, but still need to specify row numbers

ha = rowAnnotation(foo = anno_mark(at = c(1:4, 11, 28, 44,:47),
labels = genes_to_show))

## Plot heatmap
Heatmap(m, name = "mat", cluster_rows = FALSE, right_annotation = ha,
row_names_side = "left", row_names_gp = gpar(fontsize = 4))

![](https://i.imgur.com/424cPpz.png)

Best Regards,

Mohammed

Answer: Use a combination of `which()` and `%in%` to isolate the row indices that match the gene character vector:

```
mark_at = which(rownames(m) %in% genes_to_show)

ha = rowAnnotation(foo = anno_mark(at = mark_at, labels = genes_to_show))
```
Question: Hi,
I have a text file of 7500 genomic positions with their chr, start, end and want to get their nucleotide sequences. Can someone point me some tool or any thoughts on how to do it?

Thank you,

Answer: Use getfasta from bedtools : http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html
Question: Hello!

I'm sorry if this question comes of as naive or ignorant because I'm very new to Bioinformatics. I'm trying to do an alignment with STAR and was wondering if I could access a pre-made STAR index for the mm10 genome. I was told I could do this from UCSC but have had no luck finding it there.

So my question is, Are there pre-made STAR index files for the mm10 genome that I could download? And if so where and how?

Thanks in advance for any help and I'm sorry to ask such a trivial question! Let me know if there's anymore detail I can give!

Answer: I'd suggest generating your own index using the mm10 genome as per the instructions below, and using the latest gencode mouse genes. To keep things consistent (major problem in bioinformatics!!!) I'd download BOTH the genome and the annotation gtf from here [http://www.gencodegenes.org/mouse_releases/current.html][1]

You want the Comprehensive gene annotation - PRI gtf and the Genome sequence, primary assembly (GRCm38) - PRI fasta sequence (this is your genome).

[1]: http://www.gencodegenes.org/mouse_releases/current.html
Question: Given a SummarizedExperiment container, what is the quickest way to identify a batch effect from one of the covariates found within the DataFrame in the colData slot? Right now I am plotting the principal components and colouring the samples by each of the covariates. I then have to check the first few components for any separation and colour by the covariates to see which is responsible. I have a large number of libraries I have to check and was wondering if there was a Bioconductor package to perform this step? I've looked at svaseq and RUVseq but I can't see that they produce any QC plots which will tell me if an effect is present and which covariate is responsible?

Answer: Two methods are best used for analysis of batch effects.

1. PCA with annotation - as you are doing, but relies on manual visual analysis
2. PVCA - https://www.bioconductor.org/packages/release/bioc/html/pvca.html - this fits a mixed effects model to the principle components then looks at effects of various co variates, quantitatively. It is called principle variance components analysis.

Update: I highly recommend https://github.com/dswatson/bioplotr/blob/master/R/plot_drivers.R, this function it makes a great plot for examining for batch effects and more.
Question: I've been spending quite some time on following problem: I sequenced a bacterial genome using paired-end reads (SOLiD) and I have a quite good reference sequence. My goal is to detect changes in the sequenced sample compared to the reference sequence.

The detection of SNPs and small indels (a couple of bp long) is quite straightforward using the standard tools (SAMtools, GATK). However I'm stuck on the task of detecting larger Indels (tens to hundreds of bp). I tried several software and stuck with Pindel (upon a recommondation on this forum).

Because I didn't know wether to trust Pindels output, I started to simulate some data (introducing indels of several sizes into the reference), mapping the original data to the reference and checking wether Pindel was able to detect the changes. Pindel is very sensitive and could detect most of those indels, however its sensitivity is also the main problem. I find it quite impossible to differentiate between true indels and false positives. There is no good statistic regarding the significance of an observation other than the raw number of supporting reads.

My questions:

What does someone having more experience in this kind of work recommend? Any other software tools? Another approach? Or will I have to accept the fact that paired end short reads are not optimal to answer this kind of question?

For the next time: What sequencing approach would you recommend? 454 reads with de novo assembly and subsequent comparison of the contigs to the reference? PacBio? Does Illumina offer a better approach?

Thanks for any help!

Answer: Read this review paper from 2011 on Structural variant calling:

http://www.ncbi.nlm.nih.gov/pubmed/21358748

Basicly there are 3 signals you can use for structural variant calling:

1) discordant pair signal

2) readdepth signal

3) split read signal or ( with as a special case denovo assembly split contig mapping signal)

The discordant pair signal and readdepth signal you can get from paired sequencing data produced on all platforms. To use the split read signal you nead long reads and do split mapping of these reads, this is not really usefull on solid data or other short read sequences.

A good discordant pair signal SV caller is breakdancer.

A good split read signal SV caller is pindel.

A good readdepth signal SV caller is cnvator.

2 upcomping multisignal SV callers are SVMiner and Lumpy
Question: Hi,

I am a newbie with scripting so I can't find an easy solution to this question by myself and I'd like to ask for some help.

I have a long list of BED files, and for each file I want to scan them row by row and if the content of a given column contains some text I am looking for (say, for example, gene name "A") I want that full row to be copied into a new, separate bed file.

I'm looking fwd to hearing your suggestions,

thanks in advance!

Answer: If Input file contains gene name in fourth column, then it will print only that line

cat input.bed | awk 'BEGIN {OFS="\t"} { if ($4 == "geneName") { print $0 }}' > outputFile

If you have many bed files and want to loop

```
for name in $(ls *bed)
do
cat $name | awk 'BEGIN {OFS="\t"} { if ($4 == "geneName") { print $0 }}' >> outputFile
done
```
Question: I have gene symbol list like following.

2-Sep
Hist2h4
5-Mar
Nupr1
7-Sep
5-Mar
15-Sep
11-Sep
Egam-1c

I have not ever seen "2-Sep" and "5-Mar" .
These gene symbol are from where? you know?

Furthermore, do you know some ID conversion application able to change the above gene symbol to Ensemble ID.

Answer: Congrats, someone used Excel on the dataset!

https://nsaunders.wordpress.com/2012/10/22/gene-name-errors-and-excel-lessons-not-learned/
Question: Hi dears,

Can someone tell me a link to plant version of biomart please? I went through this link but it seems inactive http://www.gramene.org/biomart/martview/ea115ee510b63cb6b1680c47ebec3440, in this link http://plants.ensembl.org/biomart/martview/767d21e3edcb8fac82b9d0be5410324b there is not my interest species (citrus). I need plant biomart for id conversion and getting bulk fasta

For my purpose which part I should select from pull-down menu?

http://www.biomart.org/biomart/martview/f6673304b8ab70c4db5ce6b6418f4a96

Answer: This might have something of interest: http://phytozome.jgi.doe.gov/pz/portal.html
Question: Hi,

The below is my Dmrcate run and I do get the DMRS. but how to get the CPGids that are part of these DMRS

corfit <- duplicateCorrelation(myMs, design=designMatrix, block = biolrep)
dmrCate.<-cpg.annotate("array", myMs, what="M", arraytype = "EPIC", analysis.type="differential", design=designMatrix, coef=2, block=biolrep, correlation=corfit$consensus.correlation,fdr=0.01)
dmrcoutput_all <- dmrcate(dmrCate,lambda=500, C=5, pcutoff="fdr",min.cpgs=1,mc.cores=1)
wgbs.ranges_all.snp7 <- extractRanges(dmrcoutput_all, genome = "hg19")
write.csv (wgbs.ranges_all.snp7,file="DMRCate_DifferentailRange.txt")
write.csv(dmrcoutput_all.snp7$results,file="DMRCate_ALL.txt")

how to get the 17 Cpgids from the result shown below

"","seqnames","start","end","width","strand","no.cpgs","minfdr","Stouffer","maxbetafc","meanbetafc","overlapping.promoters"
"1","chr6",31539539,31540750,1212,"*",17,2.06654851019045e-120,2.16758469969414e-149,-0.162022961391297,-0.106211449278892,"LTA-001, LTA-002, LTA-004, LTA-003"

I tried with gRanges but unable to understand how to perform it

P.S: I have posted the same in Bioconductor Forum. sorry for redundancy. Any help soon is appriciated

Answer: gr <- granges(OBJECT)
locs.ranges<-names(gr[seqnames(gr) == "chr6" & start(gr)>=31539539 & end(gr) <=31540750 ])

The above two lines I was able to subset the ids belonging to 1 location of interest. How can I loop this multiple locations present in a dataframe?
Question: Hello,

I am very new to the world of sequencing and would really appreciate your knowledge. I am trying to study a genomic region containing 5 very homologous genes and have obtained the FastQ files generated using MinION. I used the tool NanoPlot to produce a QC report but am struggling in understanding it. Particularly regarding the quality scores and quality cut-offs; I appreciate a quality threshold depends upon my application but would it be possible to explain very simply what the number of the 'mean read quality' actually represents? If I go on to filter my reads, is there a standard cut-off for ONT data?

Any information around this would be extremely helpful!

Answer: TLDR: if you can align the reads (i.e. if you have a reference genome) then you might want to filter on mapping quality, and not on read accuracy.

In NanoPlot, the mean read quality is the mean of the base call quality scores. To be entirely correct, those (Phred scale) base call quality scores are first converted to their probability of being correct, averaged, and then back to the Phred scale. A bit more about that can be found in [this blog post][1].

ONT will by default filter on a minimal quality score of 7, but that's quite arbitrary. I don't know why they went with that score. For my applications, which is structural variant detection, I don't filter at all. If I get a low-quality read, which does map reliably, and identifies a variant then I'm happy. The quality scores [match quite well with the percent identity][2] if your data is recent at least. So as such, you can estimate the error rate of a read, based on the quality score. A quality score of 7 corresponds with a ~80% accurate reads, which is not amazing. See also the image below for how the Phred score corresponds to the probability of error or the accuracy:

![enter image description here][3]

For your application I would let the aligner judge reads: if it aligns with a high mapping quality to one of those genes, and not to the others, then that's a good thing right? I don't know how similar the genes are, and how long your reads, though, and also not what your aim is.

[1]: https://gigabaseorgigabyte.wordpress.com/2017/06/26/averaging-basecall-quality-scores-the-right-way/
[2]: https://gigabaseorgigabyte.wordpress.com/2018/08/30/update-on-oxford-nanopore-basecall-quality-scores/
[3]: https://us.v-cdn.net/5019796/uploads/FileUpload/78/663145e9df43db3efe5df4d0b88cf4.png
Question: Hello

I am trying to do some statistical analysis on RNA-Seq data. I have 5 different treatments with biological duplicates for each conditions. I am trying to analyze these data for differential gene expression among different treatments. Since pairwise testing of the data in initial stages would lead to greater amount of false discovery in genes, what is the best approach to do statistical analysis on multiple groups? Would performing ANOVA analysis be a good method to start with? And what bioinformatics package could I use to perform this analysis (edgeR or BaySeq)?

Thanks

Answer: EdgeR would work for you (it's not the only one, but it's one I've used successfully for these types of experiments). It can handle more complex experimental designs in a better way than just doing a large collection of pairwise tests.

In the manual, check out the sections "More complex experiments (GLM functionality)" and "An ANOVA-like test for any differences". In the GLM mode, you can give it an arbitrary experimental design with multiple treatments, replicates, and batches, and it will do the "right thing".
Question: I am working on wheat genome, I want to do analyze comparative genome analysis of 3 varieties of wheat .I have sequence files from illumina 1.9 in fastq format. I checked the quality of reads by fastqc tool.GC content are not in normal range.(47- 49). What is the normal value of % GC for RNA seq reads ? The other question is that kmers are also not in correct range.How can I correct it.
For trimming,adaptor seq file is required,but i don't have this file..Is that possible to remove these two error?If yes then how can I do?Can I skip the trimming step and go to next step of mapping?

In this file all parameter values are correct except kmer s and GC content.Is there any need to trim it?If yes then how can I do?

file:///home/comsats-ra/fatimamphilldata/G1_cleaned_R1_fastqc.html#M11

Answer: You can also look into **DNApi : De novo adapter prediction algorithm for small RNA sequencing data**

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

de novo adapter prediction (iterative) algorithm for small RNA sequencing data. DNApi requires Python (2 or 3) under a Linux/Unix environment. DNApi accept (un)compressed FASTQ files or redirected standard input (stdin) as an input. You can simply run:

$ python dnapi.py

or

> $ | python dnapi.py -

To see the detailed usage, type:

$ python dnapi.py [-h | --help]

DNApi can predict most 3′ adapters correctly with the default parameters. However, if you want to tweak the parameters or want to run other prediction modes, see [prediction modes and parameters] (https://github.com/jnktsj/DNApi#prediction-modes-and-parameters) for more detail.
Question: Hey ya'll,

This is the error I am getting:

> scfp3 <- SCTransform(scfp3, method = "glmGamPoi", vars.to.regress = "percent.mt", verbose = FALSE)

Error in density.default(x = genes_log_gmean_step1, bw = "nrd", adjust = 1) :
'x' contains missing values
In addition: Warning message:
In vst(method = "glmGamPoi", umi = new("dgCMatrix", i = c(0L, 2L, :
NaNs produced

I think this error might be because there are zero integers in the expression matrix. For example the data frame looks like this ![enter image description here][1]

[1]: /media/images/c9e3cf8a-d999-4fce-8ef6-fe5c78a6

Could someone help me solve this, please? I have been messing around with changing Na to zeros and vice versa, no solution yet. I will keep trying.

Answer: Turns out the data was already normalized (log-transformed).

This is why I wasn't able to run SCTransform on it.

[https://github.com/satijalab/seurat/issues/4546#issuecomment-850074482][1]

[1]: https://github.com/satijalab/seurat/issues/4546#issuecomment-850074482
Question: I am making my first UCSC assembly hub to display a non-UCSC annotated genome within the browser. All is well except that I cannot work out how to add the gene annotation, which is currently in GFF3 format. I am aware that track hubs only except the "big" file versions, so presumably a bigBed version of the annotation is required. Does anyone know of a handy method of converting GFF3 to BED/bigBED? I think BED12 is required I to retain the differentiation between CDS, UTR and introns...

Thank you!

P.S. I have Googled this! https://www.biostars.org/p/85869/ is a help, but I would be interested to know if there have been developments since then.

EDIT: GTF or GFF2 can be used for gene annotation!

Answer: In the end I contacted UCSC browser directly. I got a helpful and detailed reply that I have edited to make it clearer how the necessary programs can be obtained. This is run in 64bit Linux. IMPORTANT NOTE: my question specified GFF3 as the starting format for the annotation, but it appeared to be much easy using GTF / GFF2.

**Fetch the programs**

```
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitInfo
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/gtfToGenePred
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/extractGtf.pl
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/genePredToBed
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/genePredCheck
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ixIxx
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedToBigBed
chmod +x genePredToBed genePredToBed genePredCheck bedToBigBed faToTwoBit twoBitInfo ixIxx
```

**Download Perl scripts from their GIT repository**

http://genome-source.cse.ucsc.edu/gitweb/?p=kent.git;a=tree;f=src/hg/utils/automation

```
extractGtf.pl
ensemblInfo.pl
```

**Method**

```
# Create twoBit version of genome
faToTwoBit genome.fa genome.2bit

# Get chromosome length from twoBit genome
twoBitInfo genome.2bit stdout | sort -k2rn > genome.chrom.sizes

# Convert GTF annotation to genePred format
gtfToGenePred -infoOut=infoOut.txt -genePredExt genome.gtf genome.gp

# Check the genePred output is valid
genePredCheck genome.gp

# Convert genePred format to BED format
genePredToBed genome.gp stdout | sort -k1,1 -k2,2n > genome.bed

# Convert BED to bigBed
# extraIndex required for position/search
bedToBigBed -type=bed12 -extraIndex=name genome.bed genome.chrom.sizes genome.bb

# Required for indexing step
grep -v "^#" infoOut.txt | awk '{printf "%s\t%s,%s,%s,%s,%s\n", $1,$2,$3,$8,$9,$10}' > genome.nameIndex.txt

# Create index for position/search function in browser
ixIxx genome.nameIndex.txt genome.nameIndex.ix genome.nameIndex.ixx
```