cannin/biostars_qa huggingface dataset:
Question: Hi,

Hoping someone can help with what may seem like a simple question. I would like to use ggplot2 to plot a graph with a black border but a custom fill.

Using the following code I have managed to puoulate the graph as I would like it:

# creating color palette
>cols <- c("red" = "red", "orange" = "orange", "nonsignificant" = "darkgrey", "Increased" = "#00B2FF", "Decreased" = "#00B2FF")

# Make a basic ggplot2 object
>vol <- ggplot(data, aes(x = lfc, y = pval, color = color, labels=gene))

# inserting mnaual colors as per color pallette and more
>vol +
scale_colour_manual(values = cols) +
ggtitle(label = "Volcano Plot", subtitle = "Colored by fold-change direction") +
geom_point(size = 2.5, alpha = 1, na.rm = T) +
theme_bw(base_size = 14) +
theme(legend.position = "right") +
xlab(expression(log[2]("VitD" / "Carrier"))) +
ylab(expression(-log[10]("FDR"))) +
geom_hline(yintercept = 1, colour="#990000", linetype="dashed") + geom_vline(xintercept = 0.586, colour="#990000", linetype="dashed") + geom_vline(xintercept = -0.586, colour="#990000", linetype="dashed") +
scale_y_continuous(trans = "log1p")

Gives me the plot:

![enter image description here][1]

I understand that I need to define points as shape 21 before defining both color and fill t get different colors for border and fill. So I tried adding ,shape = 21, colour = "black" into my geom_point() agreement, but this turns everything black

> vol +
ggtitle(label = "Volcano Plot", subtitle = "Colored by fold-change direction") +
geom_point(size = 2.5, alpha = 1, na.rm = T, shape = 21, colour = "black") +
scale_colour_manual(values = cols) +
theme_bw(base_size = 14) +
theme(legend.position = "right") +
xlab(expression(log[2]("VitD" / "Carrier"))) +
ylab(expression(-log[10]("FDR"))) +
geom_hline(yintercept = 1, colour="#990000", linetype="dashed") + geom_vline(xintercept = 0.586, colour="#990000", linetype="dashed") + geom_vline(xintercept = -0.586, colour="#990000", linetype="dashed") +
scale_y_continuous(trans = "log1p")

![enter image description here][2]

THANKS!!

[1]: https://i.imgur.com/14wKc8S.png
[2]: https://i.imgur.com/xW5STjv.png

Answer: Hi,

I assume (haven't tested it on your data), that you should merely replace *color = color* with *fill = color* in your assignment of **vol**
Question: Hi guys,

I'm totally new in the field of ChipSeq histon acetylation data analysis.

I have to analyse some ChipSeq data. As from many tutorials and on-line documentation I started generating .fastq.gz from raw sequenced data. Then the trimming with Trimmomatic and then the alignment using BWA align.

After bwa align I ended up with .sai files each for each lane, i.e. each for each .fastq.gz file.

The situation is this:

X1_10_S8_L001_R1_001.aligned.sai
X1_10_S8_L001_R1_001.fastq.gz
X2_11_S8_L001_R1_001.aligned.sai
X2_11_S8_L001_R1_001.fastq.gz
X3_13_S8_L001_R1_001.aligned.sai
X3_13_S8_L001_R1_001.fastq.gz

I have thousands of files like this.

I would like to create .sam files to finally generate .bam files.

Is there a way to loop over all the "paired" files to generate the .sam files and then the bam files?

Thank you in advance for your help

B.

Answer: If your system supports parallel the try the following command with some modifications in the file name pattern.

parallel --verbose -j 20 'bwa mem -t 10 bwa_index {1}_1.fastq.gz {1}_2.fastq.gz|samtools sort -@5 -o {1}_sorted.bam -' ::: $(cat files.txt)

Change the number after `-j` to increase or decrease the number of files to be processed per batch.

`files.txt` contains all the file names.

Without parallel.

#!/bin/bash
for file in `cat files.txt`
do
bwa mem -t 10 bwa_index ${file}_1.fastq.gz ${file}_2.fastq.gz|samtools sort -@5 -o ${file}_sorted.bam -
done

*BWA Version: 0.7.17-r1194-dirty
Question: I have both short (illumina) and (pacbio) long reads of the same sample and I want to align it to a reference genome. I know there are many aligners, but I was wondering if there is one that uses a hybrid method, i.e. uses both short and long reads to improve the accuracy. I know that unicycler, for example, does a hybrid de novo assembly. However, I am looking for a program that does hybrid alignment to a reference genome. My google search did not find anything. Any comment would be appreciated.

Answer: No such hybrid aligner exists since aligners perform read alignment independently for each fragment. That is: per read for long read, and per read pair for paired-end short read sequencing. A 'hybrid' approach would give the same results as aligning the short and long reads independently with the same aligner.

Short and long reads can be combined by 'hybrid' methods either prior to alignment (using short reads to correct long reads, or perform de novo assembly), or after alignment (reference-guided assembly/alignment adjustment/variant calling based on the aligned short+long reads).
Question: Hi

In relation to the data filtering GBS, Elshire et al. 2011 describes:

> Analyses of the 86 bp sequencing reads were based upon the unfiltered qseq files, since the filtering process that produces fastq files sometimes discarded good reads that aligned perfectly to the reference genome for at least 64 bases. Starting with the qseq files from a flow cell, we first filtered for reads that (1) perfectly matched one of the barcodes and the expected four-base remnant of the ApeKI cut site (CWGC), (2) were not adapter/adapter dimers, and (3) contained no ""Ns"" in their first 72 bases.

What does it mean "Ns"?

Answer: "N" is an [IUPAC code][1] representing any base (A,C,G,T). In sequencing, this base can be thought of as being unknown.

[1]: http://www.bioinformatics.org/sms2/iupac.html
Question: Good morning, actually I'm working for a GWAS for the metabolic syndrome, I've `.bed` `.bim` and `.fam` files also `.ped` and `.map`, I'm looking for a tool that allow me to obtain `.geno` file , because I need to do some research with EIGENSOFT and it needs some file that I can't have right now. There is some tool that is capable to convert some of my file into `.geno` file.

I need also .ind and .snp and .pca files how can i obtain these files also?

(I apologize for my English I'm not British or American)

Answer: You can still work with the PLINK binary files:

cp study.bim study.pedsnp ## Renaming so that Eigensoft can recognize file type
cp study.fam study.pedind ## Renaming so that Eigensoft can recognize file type

smartpca.perl -i study.bed -a study.pedsnp -b study.pedind -o study.pca -p study.plot -e study.eval -m 0 -l study.log

rm study.pedsnp ## removing temp files
rm study.pedsnp ## removing temp files
Question: Hello everyone!

This is my first post on here. Hope you guys can give a helping hand to me.

This semester I'm getting a computational biology course. I've been enjoying though being a biochemistry based student makes dealing with calculations and maths a bit challenging.

There's this homework I'm struggling with: http://web.iyte.edu.tr/~bilgekaracali/EE550/Homework%201.pdf

Even though I'm reading the text book I can't manage to find out how I am going to create PAM and BLOSUM matrices from the aminoacid sequence given.

Can you show me a way out as I will keep searching for such examples but solved?

Thank you and sorry if I'm breaking some kind of forum rules by posting this, I just started to feel overwhelmed.

Answer: - https://web.archive.org/web/20130726133541/http://biorecipes.com/Dayhoff/code.html
- http://www.perlmonks.org/?node_id=1049813
- http://books.google.com/books?id=C88KrD4enWkC&pg=PA68&hl=en#v=onepage&q&f=false
- http://books.google.com/books?id=bvY21DGa1OwC&pg=PA96&source=gbs_toc_r&cad=2#v=onepage&q&f=false
- https://www.biostars.org/p/81517/
Question: Hi,

I have a couple of contigs, and I would like to answer questions of the following sort

a.) Are they unique?

b.) Which organism do they come from?

c.) Do they come from a protein encoding gene?

d.) What is that protein's function etc...

I am very new to Bioinformatics (even Bio for that matter) so any help would be appreciated. Even links to further readings.
Thanks.

Answer: What about starting with a simple NCBI BLAST ? http://blast.ncbi.nlm.nih.gov/Blast.cgi
using the blastx algorithm ("useful for finding similar proteins to those encoded by a nucleotide query"). You can later use the "Taxonomy reports" to get a list of the organisms found in your query.
Question: The package upsetr was installed with conda. It **works** in my **terminal**

conda activate upsetr
conda deactivate

I was trying to invoke conda in a **Makefile** , something like

SHELL=/bin/bash
test:
conda activate upsetr && conda deactivate

but it doesn't work:

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

$ conda init

Currently supported shells are:

I can reproduce this error by piping the command into bash

$ echo 'conda activate upsetr && conda deactivate' | bash

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

$ conda init

Why does it fail ?

Answer: Everything define in your .bashrc or .profile won't work when using `| bash`

If you really want to execute on a new bash instance try

```
echo 'source .bashrc && conda activate' | bash
```

Otherwise if you to run it on your current shell try
```
eval $(echo 'conda activate')
```
Question: I have gtf file containing 7 columns from which i am going extract only gene name part from column 7th.

The 7th column contains some information including `gene_id`, `gene_name` and so on that I posted one row of column 7th below:

gene_id "XLOC_000001"; transcript_id "TCONS_00000001"; exon_number "1"; gene_name "NAC001"; oId "AT1G01010.1"; nearest_ref "AT1G01010.1"; class_code "="; tss_id

I need only `gene_name` from this column, for example "NAC001", how can I extract which I need please?

Answer: Sorry I missed the part that you only posted column 7 from your file. Either use @dschika's solution below or cut the 7th column and then use the one liner below.

$ awk -F ";" '{sub(/gene_name/,""); print $4}' your_file

If you don't need the quotes around the gene name

$ awk -F ";" '{sub(/gene_name/,""); print $4}' your_file | sed 's/"//g'
Question: I was reading a paper (GREAT improves functional interpretation of cis-regulatory regions) and they refer that hypergeometric distribution is used for assessing the proximal gene enrichment in ChIP-seq data.

> In a typical analysis, one compares the total fraction of genes annotated for a given ontology term with the fraction of annotated genes picked by proximal binding events to obtain a gene-based P value for enrichment. (a) This procedure has a fundamental drawback: associating only proximal binding events (for example, under 2-5 kb from the transcription start site) typically discards over half of the observed binding events (a).

They propose that one could extend the regulatory domain and use a binomial test for distal binding sites. Is it possible to extend the regulatory domain (the way they do it in the paper) and apply a hypergeometric test on those domains?

![nbt1630.Fig1](https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fnbt.1630/MediaObjects/41587_2010_Article_BFnbt1630_Fig1_HTML.jpg)

Answer: Hi Tonja,

The rationale for the use of a binomial test is described in the paper: http://www.ncbi.nlm.nih.gov/pubmed/20436461

In the introduction, in the second paragraph, they say "...the standard approach to capturing distal events--associating each binding site with the one or two nearest genes, introduces a strong bias toward genes that are flanked by large intergenic regions' and then they further explain that this bias leads to the generation of false positive enrichment scores.

The hypergeometric test is not biased in this way for PROXIMAL regions because they do not have such wide variation in size (the authors state they are between 2-5kb). But the same cannot be said for DISTAL regions - they have huge variation in size.

So, to avoid this enrichment, the hypergeometric test is not used...Think of it this way, if you used the same procedure in Figure 1a as you did in Figure 1b (i.e., hypergeometric test both times), then genes that had huge flanking regions would pop up in your analysis way more than genes with shorter flanking regions, due to probabilistic considerations introduced by the sheer size of the region.

So, instead they define a regulatory domain **and** the number of bases that that covers (step 2 in Figure 1b), then convert to a fraction (see Results, page 495-496). This is a much better approach than using the hypergeometric test because it is free from the type of bias they describe.

In answer to your question, yes, you can use a hypergeometric test any time you wish to test for enrichment of items. However, in this case, applying it is likely to lead to a biased (systematically inaccurately estimated) test statistic.

If you still have questions after reading this, reading the paper plus references 12,15, and 16 should clarify the issue beyond and doubt.
Question: I have performed GATK's Base Quality Score Recalibration (BQSR) on my BAM files.

The BQSR has made the per base sequence quality look very bad on FastQC. Is this normal? Should I use the BAM files that have been through BQSR and therefore have worse FastQC reports?

The FastQC per base sequence quality goes from being in the green zone ![https://ibb.co/hyv5DQy][1]

to being all in the lower red zone

![https://ibb.co/ckN1fFw][2]

The next step I will do is SNP calling. This is genomic data.

(The BAM files contain only the data that is paired, matches uniquely, and has the duplicates marked by Picard). The FastQC images are from just one of my samples but all samples have the same pattern of being all green or almost all in the green zone, to being entirely in the red zone after BQSR.

[1]: https://i.ibb.co/g3C1zq3/Screenshot-2019-06-12-at-10-19-19.png
[2]: https://i.ibb.co/hLKyQWc/Screenshot-2019-06-12-at-10-19-30.png

Answer: Apparently it is because I only have data for 1 gene and GATK says "A critical determinant of the quality of the recalibration is the number of observed bases and mismatches in each bin. This procedure will not work well on a small number of aligned reads. We usually expect to see more than 100M bases per read group; as a rule of thumb, larger numbers will work better" https://software.broadinstitute.org/gatk/documentation/article?id=11081
Question: There is some weirds situation when I using simple blast search.
Version of blast is 2.9.0+

When I use the `evalue` parameter of blastp, I think it should be a simple cut off value as the same as others.
But I found out that, even with the same query and same subject. The evalue could not extract correct sequences at all.

With 1e-38 evalue, I got these sequence without the lower one (1.60e-83)
But if I use 1e-35 evalue, I got this.

This makes me question the cut off/filtration of the blast.

The situation is also described in the following picture.
https://gist.github.com/444thLiao/5d071703f65935c71906a19d39f0187f

If there any need to download the file/db I used. please feel free to download it at [download data](https://drive.google.com/open?id=1PRNs6W2DhbQqKp7jgW7HsI09qvAKzUVn)

Answer: OK, this is a tricky issue you mention here, though one that has been extensively discussed on biostars (and other fora / blog posts )

Can't find any links to point to right now, but will look for them.

long story short: the filtering is not done at the end (== reporting stage of blast) as you might intuitively expect but also has it consequences earlier on in the blast process causing differences in output.

EDIT: here is some background info on it :

biostars posts:

- https://www.biostars.org/p/52787/
- https://www.biostars.org/p/340129/
- https://www.biostars.org/p/341227/

a blog post (which kinda started it all):

https://blastedbio.blogspot.com/2015/12/blast-max-target-sequences-bug.html

and a paper on this:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6662297/
Question: Take one VCF file for example,

##fileformat=VCFv4.1
##samtoolsVersion=0.1.16 (r963:234)
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the site allele frequency of the first ALT allele">
....
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SLVi33.16 SLVi33.25 SLVi33.26
chr1 85691 . C T 6.26 . DP=2;AF1=1;CI95=0.3333,1;DP4=0,0,2,0;MQ=40;FQ=-28.9 GT:PL:GQ 0/1:0,0,0:3 1/1:20,3,0:4 1/1:18,3,0:4
chr1 88338 . G A 6.95 . DP=3;AF1=1;CI95=0.1667,1;DP4=0,0,2,1;MQ=30;FQ=-28.1 GT:PL:GQ 0/1:0,0,0:3 0/1:0,0,0:3 1/1:38,9,0:8

The DP=2 at chr1:85691, however, there are three individuals in the file, why does every individual has genotype? As AF1 means allele frequency of ALT allele, it supposed to be AF1=5/6 at chr1:85691, why it's truely AF1=1? I'm so confused, anybody can help me? Thanks in advance.

Answer: DP is just the raw read depth at that position. DP of 2 just means 2 reads out of your 3 individual samples mapped to that position and passes the quality filter. It doesn't say whether it is 2 reads from one of your individual or 2 read from 2 individuals..etc.

AF1 of 1 pretty much just means the 2 reads (2/2) both have a T instead of an C at that position.

Why are you expecting it to be 5/6? Are you trying to compare this to some other data?

edit **

Samtools is really not great at outputting and calculating the genotype field. Here is another biostar questions about it.

The PL fields showing 0,0,0 pretty much means there are no reads covering it. It looks like samtools falsely calculates 0,0,0 as heterozygous because the probabilities are similar (all zeroes in this case).
Question: Is there a tool like bedtools shuffle which I can use to randomly shuffle a bam file? or will I have to convert my bam into a bed and then shuffle it? thanks

Answer: `samtools bamshuf file.bam` does the job. It'll also break up your bam in `-n` (default=64) chunks. Unfortunately though, it wants to shuffle the entire bam and then writes the chunks to disk (does not stream even with the `-O` option). If you just want to sample a few reads from a large bam quickly, you can "shuffle" some coordinates and then use `samtools view file.bam chr:start-start | head -1` to get a read. If you want read pairs, you can either use `samtools bamshuf` or sort the bam by name(`samtools sort -n`).
Question: Hello, sorry if this is a really stupid question but I have been struggling with this for a week or so now.

Basically I carried out a bacterial RNA-seq on 2 conditions in triplicates. After doing some analysis using DESeq2 it is looking like there is no differential expression between the 2 conditions however, is there a way that I can save the normalised count data to a cvs file so that I can check the quality of the replicates?

Thanks

Answer: foo <- counts(object, normalized = TRUE)
write.csv(foo, file="norm_counts.csv")

where `object` is a `DESeqDataSet` object
Question: Hi,

I was working on converting my vcf file into plink ped/map format.

I used the following command from vcftools

vcftools --vcf myvcf.vcf --plink --out myplink

I also would like to filter my snps based on minor allele frequency (maf).

I was wondering whether simply using

vcftools --vcf myvcf.vcf --plink --maf 0.05 --out myplink

would do it in one step?

Or shall i do it in two steps. First convert it to the ped/map format and then use plink to do the filtering?

Answer: PLINK 1.9 can do this:

plink --vcf myvcf.vcf --maf 0.05 --recode --out myplink

(Are you sure you want to convert to ped/map, though? PLINK prefers to work with binary bed/bim/fam.)
Question: Hi.

I know that what soft-clipped reads are.

But when facing sam files, also, there are many hard clipped reads.

and I cannnot figure out what differences between them even though I check the sam specification.

Could you explain the difference of them?

I attached the easily understood figure

![< image not found >][1]

[1]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2723002/figure/F1/

Answer: soft-clipped: bases in 5' and 3' of the read are NOT part of the alignment.

hard-clipped: bases in 5' and 3' of the read are NOT part of the alignment AND those bases have been removed from the read sequence in the BAM file. The 'real' sequence length would be length(SEQ)+ count-of-hard-clipped-bases
Question: Hello,

I was given the abundance files from kallisto for a project. However, I also needed the genomebam files which where given to me later. I noticed that for the bam files the kallisto command was run differently. The first files were run without `-b` but when the bam files were generated they were run with `-b 100`. Would the bootstrap option change anything on the BAM files? I just want to visualize them on IGV to check if my knock out has worked. I am not sure if I can use the BAM files I was provided with or if I should ask for them again.

Thank you

Answer: No, bootstraps have no effect on BAM files.
Question: I want to use a binary matrix to build a phylogenetic tree using the R package [phangorn][1] as done in this paper:[Tracking the Genomic Evolution of Esophageal Adenocarcinoma through Neoadjuvant Chemotherapy][1]

In the methods they say:

> Trees were built using binary presence/absence matrices built from the regional distribution of variants within the tumor. The R Bioconductor package phangorn (1.99-7; ref. 36) was utilized to perform the parsimony ratchet method (18), generating unrooted trees. Branch lengths were determined using the `acctran` function.

I have the binary presence/absence matrix, however the phangorn package uses phyDat objects, which are derived from sequence alignments according to the [phangorn vignettes][3].

My question is:

**How can I use a binary table to build a phylogenetic tree with the R phangorn package?**

If there is a way to read the binary matrix as a phyDat object, that would solve the problem, but I don't see how that could be done.

[1]: https://github.com/KlausVigo/phangorn
[2]: http://cancerdiscovery.aacrjournals.org/content/5/8/821.long
[3]: https://cran.r-project.org/web/packages/phangorn/vignettes/Trees.pdf

Answer: Hello Alejandro,

There are generic functions `as.phyDat()` in phangorn to transform matrices and data.frames into phyDat objects.

For example you can read in your data with `read.table()` or `read.csv()`, but you might need to transpose your data. For matrices `as.phyDat()` assumes that the entries each row belongs to one individual (taxa), but for data.frame each column. For binary data you can transform these with a command like (depending how you coded them):

as.phyDat(data, type="USER", levels = c(0, 1))
as.phyDat(data, type="USER", levels = c(TRUE, FALSE))

Regards,
Klaus