cannin/biostars_qa huggingface dataset:
Question: I have a data frame with the differentially expressed genes from EdgeR, Now I am trying to make a volcano plot of it but I want to see only selected genes that are of interest to me to be labelled on the volcano plot.
My data frame looks like this

head(results)

Gene Fold pvalue FDR sig
1 ADORA2A 10.273854 1.164636e-28 2.234471e-24 FDR<0.05
2 IL23A 8.132293 4.554177e-28 4.368822e-24 FDR<0.05
3 IL1A 8.430078 6.768343e-27 4.328581e-23 FDR<0.05
4 CXCL6 7.102900 3.299464e-23 1.582588e-19 FDR<0.05
5 CCR7 8.950486 9.111421e-23 3.496235e-19 FDR<0.05
6 IL18R1 6.759646 7.283440e-22 2.329001e-18 FDR<0.05

I tried the below code to generate the volcano plot, which is generated successfully and the selected genes are marked too but not with their exact name but with some random number.

results$genelabels <- ""
results$genelabels <- ifelse(results$Gene == "IL23A"
| results$Gene == "IL1A"
|results$Gene == "IL6"
|results$Gene == "CD80"
|results$Gene == "CD86"
|results$Gene == "NFKB"
|results$Gene == "BAFT2", TRUE,FALSE)
ggplot(results) + geom_point(aes(Fold, -log10(FDR),col=sig))+ geom_text_repel(aes(Fold, -log10(FDR)),label = ifelse(results$genelabels == TRUE, results$Gene,""), box.padding = unit(0.45, "lines"),hjust=1) + theme(legend.title=element_blank(),text = element_text(size=20))+ scale_color_manual(values=c("red", "black"))

Can somebody help me out with what is wrong with the code that is does not show the exact name of the gene on volcano plot?

Answer: Your labels are of class **factor** and not **character.** If you wrap it with `as.character` in your plotting call it should work.

ggplot(results) + geom_point(aes(Fold, -log10(FDR),col=sig))+ geom_text_repel(aes(Fold, -log10(FDR)),label = ifelse(results$genelabels == TRUE, as.character(results$Gene),""), box.padding = unit(0.45, "lines"),hjust=1) + theme(legend.title=element_blank(),text = element_text(size=20))+ scale_color_manual(values=c("red", "black"))

Specifically change the following:

ifelse(results$genelabels, results$Gene,"")

to

ifelse(results$genelabels, as.character(results$Gene),"")
Question: Hello,

How can I ignore the first bases of index reads with bcl2fastq2 `--use-bases-mask` option?

The index read has 8 cycles, but I want to use only "last 7" of them. illumina bcl2fastq manual only mention about ignoring the last reads, not ignoring the first reads.

I tried the following option but it didn't work.

--use-bases-mask y151,nI7

Answer: You were very close:

--use-bases-mask y151n,I7
Question: Hi - I am looking over some raw data from a recent RNA-seq run, and when I put the FASTQ files through FastQC, everything looked normal (i.e. normal output from the analyses) except the per-base sequence quality. This is the graph I got for one of my 8 samples - all of them looked similar (not identical), though:

![weird sequence quality plot][1]

This is my first time analyzing my own raw data from our core facility, so of course it is possible that they are just so good that they have essentially no read-to-read variation in sequence quality. But I would have expected something more like this (exact same analysis, on a FASTQ file downloaded from SRA):

![enter image description here][2]

Before I talk to the core people, I was wondering if anyone here had ever encountered anything similar, and if so, what could cause it. Thanks a lot!

[1]: http://imgur.com/i24hdmN.png
[2]: http://i.imgur.com/aIReovC.png

Answer: We have seen very similar output with the newer chemistries for RNA-Seq data. It's normal ![enter image description here][1]

[1]: http://i.imgur.com/RQJ0wB1.png
Question: Hi,

I have a raw counts file extracted by featurecounts tool and also a gene length matrix. I was going to count rpkm for my 76 samples then I did like below

```r
rawcounts <- read.table("rawcounts.txt", header = T, sep = "\t")
head(rawcounts[,1:4])
Geneid SRR018347 SSR019035 SRR074122
1 AT1G01010 31 37 4
2 AT1G01020 153 211 19
3 AT1G01030 10 13 0
4 AT1G01040 88 106 3
5 AT1G01046 0 0 0
6 AT1G01050 228 272 195

dim(rawcounts)
[1] 33602 77

gene-length <- read.table("gene-length.txt")
head(gene-length[,1:4])
V1 V2 V3 V4
V1 Length 1688 1774 1905

dim(gene-length )
[1] 1 33603

d <- DGEList(counts=rawcounts)
cpm(d)
cpm(d,log=TRUE)

d$genes$Length <- gene-length
rpkm(d)

Warning message:
In Ops.factor(left, right) : '/' not meaningful for factors
```

Even I tried with transposed gene length matrix but again I got error

What is my fault please?

Answer: The short code I used to counts rpkm with edgeR package

```r
# rpkm normalization by edgeR package
# I need a file containing gene length that I extracted from featurecounts tool output
# reading gene length file
mycounts1 <- read.csv("length.txt", stringsAsFactors = FALSE)
# read as numeric
mycounts1 <- as.numeric(mycounts1$Length)
# transposing the gene length file
tmycounts1 <- t(mycounts1)
# setting as numeric
sapply(tmycounts1, class)
# reading the rawcounts file
mycounts <- read.table("rawcounts.txt", header = T, sep = "\t")
# watching the head of uploaded file
head(mycounts[,1:4])
# watching the dimension of matrix
dim(mycounts)
# assigning the column one as mycounts rowname
rownames(mycounts)<- mycounts[,1]
# read as matrix
as.matrix(mycounts)
# viewing mycounts class
class(mycounts)
# setting as numeric
sapply(mycounts, class)
# removing the column on
mycounts <- mycounts[,-1]
# watching the head of rawcounts file
head(mycounts[,1:4])
# applying edgeR function
library(edgeR)
# assigning the rawcounts to a new variable
d <- DGEList(counts=mycounts)
# reading gene length file as a matrix
d$genes$Length <- c(tmycounts1)
# rpkm function
mycounts <- rpkm(d)
# writing the result
write.table(mycounts, file = "rpkm.txt", dec = ".", sep = "\t", quote = FALSE)
```
Question: I have been trying to find if TCGA PRAD (prostate cancer) RNA-seq data strand specific or not. I can not find any relevant info. I did find that it is Illumina protocol, but nothing more. Does anyone if they are strand specific? Do you have the right link?
Thanks

Answer: All TCGA samples are non-stranded.
Question: How to convert .fna.gz files to .fna in a folder?

I would also like to convert all .fna to .fasta.

Thanks in advance

Answer: I hope this will work.

# Rename all *.fna to *.fasta
for f in *.fna; do
mv -- "$f" "${f%.fna}.fasta"
done

Or you can try

find . -name '*.fna' -exec rename 's/\.fna$/.fasta/' \{} \;
Question: I have multifasta sequence and i am trying to rename the header and print the first id only

>gene_1|id=abc|protein=capA|ABD5672.4
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTCG
TACGGGTTTT

desired output:

>gene_1
ATGCGTCGACGTCGTACGGGTTTT
CGTACGGGTTATGCGTCGACGTC
GTACGGGTTTT

I tried following using sed but didnt work it seems sed -r 's/^([^|]+\|){0}//;s/\|//'any improvement on this? or any other alternative ?

Answer: ```cat gene.fasta | paste - - |awk '{ split($1,a,"|");print a[1] "\n" $2}'```
Question: I have a list of genes i.e.:

FOXG1
SNRNP70
CREBBP
JUN
FHL1
SNW1
FHL3
FHL2
...
...

My question is how can I use python or perl to automatically retrieve the corresponding accession numbers for the UniProt database? I'm looking for the name of the package and classes/methods as I'll be able to figure out the programming myself. I'm just having difficulties identifying the correct tools for the job.

Thanks

Answer: You can use one of the APIs for biomart. Here's an example query, just click on the perl link for how to reproduce it with the perl API.
Question: I can't find much info on this question online, so apologies if it is answered elsewhere. I'm working with whole exome seq data with the aim of identifying somatic variants. Using fastqc I can see adaptor in around 20% of reads by the end of the read (150bp PE).

I always understood that adaptor trimming was routine in this situation but I have seen a few resources recommend against it, without providing any explanation as to why. Can anyone explain why it might be detrimental to trim adaptors prior to alignment?

Thank you

Answer: It is not detrimental to trim adapters. Reason some people choose not to do it is to save time, especially if you are aligning to a good reference (which you may be doing). Most aligners will soft-clip parts of read that do not align so in theory can take care of any adapter sequences that way. Personally I take the extra time to remove them so I no longer need to think about their presence for any downstream analysis.
Question: I am using 'bwa' to map Illumina reads to a large database and would like to be able to store (and analyse) the information on multiple hits. I have forced this very brutishly using the following:

bwa samse -n 1000 database.fasta aln_sa.sai reads.fasta > aln_1000.sam

However, the output is no longer in classical .sam format (i.e. no flags, no CIGAR string, ...). Would someone know how to force a real .sam output?

Just as an example, here is the kind of output I get:

>1-815511 3 3
bta-mir-21 +8 0
mRNA___ENSBTAG00000029897|ENSBTAT00000042276|bta-mir-21|bta-mir-21 +18 0
19 +11033079 0
>2-592417 3 3
bta-mir-21 +8 0
mRNA___ENSBTAG00000029897|ENSBTAT00000042276|bta-mir-21|bta-mir-21 +18 0
19 +11033079 0
...

What I have tried so far with no luck:

use the xa2multi.pl script which parses the information in the XA:Z: (information on other hits when their number is inferior to the value of the -n option) to output one line per hit -- this scripts is now part of the 'bwa' package.

Unfortunately, the output of a normal:

bwa samse database.fasta aln_sa.sai reads.fasta > aln.sam

does not contain any XA:Z: field... Even for read 2-592417 which has several hits:

1-815511 0 mRNA___ENSBTAG00000029897|ENSBTAT00000042276|bta-mir-21|bta-mir-21 18 0 23M * 0 0 TAGCTTATCAGACTGATGTTGAC * XT:A:R NM:i:0 X0:i:3 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:23
2-592417 0 19 11033079 0 22M * 0 0 TAGCTTATCAGACTGATGTTGA * XT:A:R NM:i:0 X0:i:3 X1:i:0 XM:i:0 XO:i:0 XG:i:0MD:Z:22

install the previous bwa-0.5.7 version for which a patch doing exactly what I need has been developed at Broad Institute. But the patch has vanished from its ftp site and is now unavailable.

Any ideas?

Answer: Thanks to Istvan's suggestion, in the latest bwa version, the output is in a fully populated SAM file (XA:Z: field is now present). I just then run the xa2multi.pl on it to obtain one hit per line.
Question: I have a gff file containing genes predicted by AUGUSTUS, the file already containing CDS, exons, and protein sequences.
I need to extract protein sequence from the file using bash.

Answer: Oh, I see. Strange, and more complicated output. Not impossible though. I'm assuming you want the hash (#) and space removed from the beginning of each protein sequence as well.

Let me know if this does the trick for you:

awk '/# protein sequence/{a=1}/# Evidence/{a=0}a' Genes.gff | sed 's/# //'

- /# protein sequence/ matches lines having this text, as well as /# Evidence/ does.
- /# protein sequence/{a=1} sets the flag when the text # protein sequence is found.
- /# Evidence/{a=0} unsets the flag when the text /# Evidence is found.
- The final a is a pattern with the default action, which is to print $0: if flag is equal 1 the line is printed.
- Finally, sed removes the hash and space from the beginning of each line.
Question: Hi everyone! :)

I have a question and did not find any answer by personal search. I would like to make a heatmap with categorical variables (a bit like this one: [heatmap-like plot, but for categorical variables][1]), and I would like to add on the left side a phylogenetic tree (like this one: [how to create a heatmap with a fixed external hierarchical cluster][2]). The ideal would be to adapt the second one since it looks much prettier! ;)

Here is my data:

A newick-formatted phylogenetic tree, with 3 species, let's say:

((1,2),3);

A data frame:

```
x<-c("species 1","species 2","species 3")
y<-c("A","A","C")
z<-c("A","B","A")
df<- data.frame(x,y,z)
```

(with A, B and C being the categorical variables, for instance in my case presence/absence/duplicated gene).

Would you know how to do it?

Many thanks in advance!

[1]: http://stackoverflow.com/questions/12998372/heatmap-like-plot-but-for-categorical-variables
[2]: http://stackoverflow.com/questions/15153202/how-to-create-a-heatmap-with-a-fixed-external-hierarchical-cluster

Answer: I figured out to do it! Here is my script for those that are interested:

```
#load packages
library("ape")
library(gplots)

#retrieve tree in newick format with 3 species
mytree <- read.tree("sometreewith3species.tre")
mytree_brlen <- compute.brlen(mytree, method="Grafen") #so that branches have all same length

#turn the phylo tree to a dendrogram object
hc <- as.hclust(mytree_brlen) #Compulsory step as as.dendrogram doesn't have a method for phylo objects.
dend <- as.dendrogram(hc)
plot(dend, horiz=TRUE) #check dendrogram face

#create a matrix with values of each category for each species
a<-mytree_brlen$tip
b<-c("gene1","gene2")
list<-list(a,b)
values<-c(1,2,1,1,3,2) #some values for the categories (1=A, 2=B, 3=C)
mat <- matrix(values,nrow=3, dimnames=list) #Some random data to plot

#plot the heatmap
heatmap.2(mat, Rowv=dend, Colv=NA, dendrogram='row',col =
colorRampPalette(c("red","green","yellow"))(3),
sepwidth=c(0.01,0.02),sepcolor="black",colsep=1:ncol(mat),rowsep=1:nrow(mat),
key=FALSE,trace="none",
cexRow=2,cexCol=2,srtCol=45,
margins=c(10,10),
main="Gene presence, absence and duplication in three species")

#legend of heatmap
par(lend=2) # square line ends for the color legend
legend("topright", # location of the legend on the heatmap plot
legend = c("gene absence", "1 copy of the gene", "2 copies"), # category labels
col = c("red", "green", "yellow"), # color key
lty= 1, # line style
lwd = 15 # line width
)
```

And I don't know how to show the result but it does work ;)
Question: Hi,

Is any tool/method to convert article title to pubmed ID? Anyone know?

Thank you,

Answer: You should be able to use [eutils][1] to query Pubmed with the title and get the PubmedID.

Also, try using [biomaRt][2] or [reutils][3]. They will help you work with the query output from NCBI as objects, which you can then use for further downstream analyses.

You'd just have to replicate this eutils query if you use reuils:

http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=pubmed&term=p.R301X%20Mutation%20and%20Variable%20Phenotypic%20Appearance%20of%20Fabry%20Disease&retMode=json

This searches for the paper titled "**p.R301X Mutation and Variable Phenotypic Appearance of Fabry Disease.**", something I picked at random.

[1]: http://www.ncbi.nlm.nih.gov/books/NBK25500/
[2]: https://bioconductor.org/packages/release/bioc/html/biomaRt.html
[3]: https://cran.r-project.org/web/packages/reutils/reutils.pdf
Question: Hi

I have a folder with 4000 files in it. I would like make a new folder, and copy into it all of the files which contain >50 fasta sequences. How do I do this?

I know that I need to create a simple loop, then grep '>' | wc -l and select only those with >50. But I am new to programming and am unsure how to write this properly.

Answer: This is a very interesting question, but you should implement it without a FOR loop, which are a bad thing.

For example you can do it with find:

mkdir -p smallfiles_folder/
find . -type f -print0 | xargs -0 grep -m 50 -o -H -c '>' | gawk -F":" '$2 > 2 {print $1} ' | grep -v 'smallfiles_folder' | xargs -i cp {} smallfiles_folder/

Explanation:

# find all files in current folder. Add -iname "*fasta" to restrict to fasta files only
find . -type f -print0 |

# count number of fasta headers in file. The -m 50 option is used to stop grepping after 50 matches, saving some calculation time) \
xargs -0 grep -m 50 -o -H -c '^>' |

# Use files to select files with less than 50 matches
gawk -F":" '$2 < 50 {print $1} ' |

# Remove target folder (avoid potential recursive loops)
grep -v 'smallfiles_folder' |

# copy / move files to target folder
xargs -i cp {} smallfiles_folder/
Question: I'm using Bowtie2 for mapping reads to human genome. I'm using Bowtie2 first time, so what are the best parameters to set if we are using these mapping results to further use to detect SNPs? bowtie2 has also given some preset options. which is more appropriate for finding SNPs from alignment results.

Which one is more appropriate local or end-to-end mode? if my purpose of alignment is SNP detection?

Answer: After what Devon say about trimming

What I am using

$BT2_HOME/bowtie2 --local -x specious -U $BT2_HOME/example/reads/longreads.fq -S eg.sam

`-U` Comma-separated list of files containing unpaired reads to be aligned, e.g. `lane1.fq,lane2.fq,lane3.fq,lane4.fq`. Reads may be a mix of different lengths. If `-` is specified, `bowtie2` gets the reads from the "standard in" or "stdin" filehandle.

`-x` The basename of the index for the reference genome. The basename is the name of any of the index files up to but not including the final `.1.bt2` / `.rev.1.bt2` / etc. `bowtie2` looks for the specified index first in the current directory, then in the directory specified in the `BOWTIE2_INDEXES` environment variable.

`-S` File to write SAM alignments to. By default, alignments are written to the "standard out" or "stdout" filehandle (i.e. the console).

If your data not indexed use

$BT2_HOME/bowtie2-build -f ~/path/to/your/genome.fa ~/where/you want/index/to/be/nameOfTheIndex

summary of the parameters: http://chipster.csc.fi/manual/bowtie2-paired-end.html

for more visit http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
Question: Hi i am beginner Bioinfo

I analysis linkage disequilibrium(ld), from plink

plink.ld result below that

CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2
3 30780810 . 3 30780827 . 0.996637
3 30780810 . 3 30780907 . 0.479487
3 30780810 . 3 30781157 . 0.715143
3 30780810 . 3 30781205 . 0.692494
3 30780827 . 3 30780907 . 0.468367

I wanna try draw LD heatmap using plink.ld result file like haploview plot

Please help me

Answer: There are quite a few old threads on this topic, but some out of date. I tried to aggregate some of these methods in my recent answer, here: https://www.biostars.org/p/379447/#379464

Essentially, you can export your PLINK data to HaploView format, and then easily compute LD and generate a LDheatmap. Given the data that you have, currently, you could probably generate the LDheatmap manually, too.

Please see the related StackOverflow thread, too: https://stackoverflow.com/questions/56118669/manhattan-plot-and-linkage-disequilibrium-heatmap-on-one-plot

Kevin
Question: I understand that pathway gene lists in databases also include the inhibiting genes.
Suppose I have a sorted list of genes, from high differentially expressed genes (hyper expression) through neutral genes to negatively expressed genes (hypo expression). When searching for pathways, should I submit one list with **both hyper and hypo subsets**? Are there insights to be gained when submitting each subset individually?

Answer: What kind of pathways are you talking about? Please define. A metabolic pathway? Kinase pathway? As you can see already with my two examples, pathways can be very different and therefore differently regulated. In metabolic pathways, you'll have kinetics with enzymes, more expression of a gene might correlate with more of that particular enzymatic reaction, however some enzymes may also inhibit a reaction (so down regulation of some members of the pathway could be expected). So I would say, both up- and down regulated genes would be interesting is such a pathway. However for a kinase pathway, that is regulated by post-translational modifications, looking at up- or down genes wouldn't really make sense. So hard to give one answer, it depends on your definition of pathway.
Question: Hello All,

I am running variant calling on some species whose reference genomes have a very high number of contigs (sometimes >400,000). The variant caller I am using splits the job by the number of chromosomes, and is overwhelmed when this number is too high. Therefore I would like to concatenate the contigs for a given species reference fasta file into ~30 contigs.

I believe I could use code such as below to merge all the contigs into one:

> grep -v "^>" test.fasta | awk 'BEGIN { ORS=""; print
> ">Sequence_name\n" } { print }' > new​.fasta

However I would like to merge them into ~30 contigs so the process can still be parallelised. I would also like to insert 1000 'N' characters between each of the merged contigs within these merged contigs, to avoid mapping issues that could be caused by merging contig sequences from different parts of the genome.

Does anyone have any advice for how to do this or know of any application that could do something similar?

Thanks in advance for your help.

Answer: You can try this `awk` solution. But first you have to calculate how many contigs should be merge into a single sequence.

Save this script as `merge_contigs.awk`

function print_n(n)
{
text = ""
for(j=0;j" {
if(n<=i) {
header = header""substr($0,2)" ";
if(n > 1) seq=seq""print_n(1000)
n++;
}
else {
print ">"header"\n"seq;
n=1;
header=substr($0,2);
seq=""
}

next;
}
{seq=seq$0}
END {print ">"header"\n"seq;}

Run it like this:

$ awk -v i=140 -f merge_contigs.awk input.fa > output.fa

Set the value for `-v i` to the number of contigs that should get merged into a single sequence.

fin swimmer
Question: Good morning.

I am trying to use kmergenie v1.6972. But it is broken when I typed "make"

```
sudo make
g++ -o minia/Pool.o -c minia/Pool.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
make: g++: Command not found
make: *** [minia/Pool.o] Error 127
```

Firstly I install g++ like the next user told me. But it appears new errors. I have r-base, Rscript installed.

I show the error with `make`, `make check` and `./kmergenie`

```
make
g++ -o minia/Pool.o -c minia/Pool.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
g++ -o minia/Bank.o -c minia/Bank.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
In file included from minia/Bank.cpp:21:0:
minia/Bank.h:13:63: fatal error: zlib.h: No such file or directory
#include // Added by Pierre Peterlongo on 02/08/2012.
^
compilation terminated.
make: *** [minia/Bank.o] Error 1
```

```
./kmergenie /home/zhu/File.fq
running histogram estimation
Execution of specialk failed: [Errno 2] No such file or directory
fitting model to histograms to estimate best k
could not predict a best k value
Execution of decide failed (return code 0). If this is a fresh Kmergenie install, try running 'make check'.
```

```
make check
g++ -o minia/Bank.o -c minia/Bank.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
In file included from minia/Bank.cpp:21:0:
minia/Bank.h:13:63: fatal error: zlib.h: No such file or directory
#include // Added by Pierre Peterlongo on 02/08/2012.
^
compilation terminated.
make: *** [minia/Bank.o] Error 1
```

Regards

Answer: Whenever you see an error such as

XXX.h: No such file or directory

the solution is to search for and install a package usually (but not always) named something like `XXX-dev`.

So in this case

```
apt-cache search zlib dev
sudo apt-get install zlib1g-dev
```

It is always worth doing a Google search with errors like these; they are common and many people will have posted solutions to them.
Question: Hi, apologies for this basic questions, I new in NGS quality control. I have been check my NGS data (Illumina - HiSeq 2500 2*100pb) using FastQC after trimming Nextera Adaptater with bbduck (BBTool package) using trimming overlap (ktrim=r k=25 mink=11 hdist=1 tpe tbo). And the checks fails Per base sequence content, Per sequence GC content and Kmer Content. My question is should I be trying remove the first base (~15) of my sequence? When I try this, Kmer Content fails again but at the end of the sequence. Should I remove again at the end?
After my goal is to assemble them with a De Bruijn assembler and thus i'm afraid of over-represent kmer.

Thanks for the help.

![enter image description here][1]
![enter image description here][2]

[1]: https://i.imgsafe.org/498873a.png
[2]: https://i.imgsafe.org/8d597c0.png

Answer: In fragmentation the Nextera kit enzyme prefers to cut DNA from certain places, which is probably the reason for the bias you observe. I would let it be..