cannin/biostars_qa huggingface dataset:
Question: I am trying to install bcl2fastq from source.

However, during the make process, I am encountering an error:

make[2]: *** [c++/lib/demultiplex/CMakeFiles/casava_demultiplex.dir/BclDemultiplexer.cpp.o] Error 1
make[1]: *** [c++/lib/demultiplex/CMakeFiles/casava_demultiplex.dir/all] Error 2
make: *** [all] Error 2

Does anyone know how you can get around this, and what this means? I have built the boost that comes with the redist folder and have pointed to this when configuring the bcl2fastq program.

Answer: Please download bcl2fastq2-v2.17.1.14-Linux-x86_64.rpm from above link and follow the instructions below: (copy/pasted from: https://github.com/igorbarinov/bcl2fastq/blob/master/install-2.17.sh)

sudo apt-get install alien --assume-yes
sudo alien bcl2fastq2-v2.17.1.14-Linux-x86_64.rpm
sudo dpkg -i bcl2fastq2_0v2.17.1.14-2_amd64.deb

I was able to generate deb (bcl2fastq2_0v2.17.1.14-2_amd64.deb) and install it. Let us know if you are not able to generate deb. bcl2fastq is available in /usr/local/bin.

my system:

$ cat /etc/lsb-release
DISTRIB_ID=LinuxMint
DISTRIB_RELEASE=18.2
DISTRIB_CODENAME=sonya
DISTRIB_DESCRIPTION="Linux Mint 18.2 Sonya
Linux mint 18.2 is based on ubuntu xenial (16.0.4)

$ uname -a
Linux genomics 4.8.0-53-generic #56~16.04.1-Ubuntu SMP Tue May 16 01:18:56 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux

For some reason, if you fail to generate deb, i can upload the deb file and share a link to download.
Question: I want to remove the decimal from the ensembl gene ID ,since it contains the decimal point it becomes difficult when i try to map the same to gene name .

gene	"nH1.bam"	"nH2.bam"	"nH3.bam"	"nH4.bam"
"ENSG00000238164.4"	-0.6534833425	-0.6404869759	-0.5898568965	-0.586357257
"ENSG00000049249.6"	1.0589150487	0.2235087421	0.5028436068	0.5201173416

I want this in my gene field "ENSG00000049249" instead of this "ENSG00000049249.6"

I tried this `awk '{gsub(/\..*$/,$1)}1'` it seems it messing up the data frame im not sure what im doing wrong.

Any help or suggestion would be highly appreciated

Answer: sed 's/\(ENSG[0-9]*\)\.[0-9]*/\1/g' input.txt
Question: I use picard software to mark duplicates, and here is my command:

```
java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms8g -Xmx16g -Djava.io.tmpdir=tmp -jar ./picard.jar MarkDuplicates I=input.bam O=out_markdup.bam METRICS_FILE=out.metrics ASO=coordinate VALIDATION_STRINGENCY=LENIENT
```

It works well but when the `input.bam` file gets bigger, the speed is very slow! I found that the picard MarkDuplicates doesn't support multiple threads. So, is there anyway to speedup picard? Another way, is there any better software to do the same as picard MarkDuplicates but with less time? I know [elprep](https://github.com/ExaScience/elprep) is another choice, but it needs very large memory!

Besides, I found that samtools can also remove duplicates, but according to my search, samtools can not remove the duplicates cross different chromosomes, so picard is more universe.

Any reply will be much appreciated!

Answer: Use `clumpify.sh` from BBMap suite (https://www.biostars.org/p/225338/ ). You do not need to align the data. Clumpify can address all types of duplicates (PCR, optical and so on).
Question: Hi all,

For plotting heatmap, I used pheatmap just by taking numerical values but not sure how to annotate and to divide heatmap based on column category.

Code which I used just taking numeral value;

x <- read.table("pheatmap.txt",sep="\t",head=T,row.names=1,check.names=F)
library('pheatmap')
pdf("16s.pdf", height=20,width=13)
pheatmap(x,scale="row",color =colorRampPalette(c('blue','yellow'))(12),cluster_cols =F,show_colnames = T,cluster_row=T,border_color="NA", fontsize_row=10,fontsize =10)
dev.off();

Now Could you please help how should I take only numerical value from the below file for pheatmap and then how to divide the plot into three different category; leaf, root, soil, and then want to add a color bar on top based on season?

Season	dry	dry	wet	wet	dry	dry	wet	wet	dry	dry	wet	wet
Compartment	Leaf	Leaf	Leaf	Leaf	Root	Root	Root	Root	Soil	Soil	Soil	Soil
Burkholderiaceae0	0	0	0.103	0.244	0.001	0.143	0.08	0.032	0.01	0.01	0.032
Bradyrhizobium	0	0	0	0	0.1354	0.094	0.12	0.16	0.211	0.005	0.00744	0.006
Anaeromyxobacter0.04	0	0	0	0.0874	0.139	0.058	0.07	0.11	0.06	0.026	0.034
Geobacter	0.001	0.0004	0	0	0.01	0.016	0.0146	0.08	0.016	0.05	0.0283	0.02
Microbacterium	0.31	0	0.02	0.038	0.003	0.03	0.0038	0.0	0	0.002	0.008	0.002
Rhizobiales	0	0	0	0	0.008	0.001	0.011	0.01	0.0121	0.01	0.4	0.0164

I would like to get image like this and also shown in this publication Fig1 (https://www.nature.com/articles/s41396-021-00974-2); instead of soil type I want to have Compartment; leaf, root, soil and then season (wet/dry) at the top

Many thanks,

[1]: https://i.ibb.co/1994gFN/GENERA.jpg

Answer: I would approach this problem in the following way: **first**, convert the column names in your count matrix to sample names; **second**, generate a second data.frame with the factors you want to divide the plot by. Then, you can use the ``annotation_col`` parameter to divide the plot in the way you want.

A simplified example of this is shown below:

> matrix <- structure(list(S1 = c(0, 0, 0.04, 0.001, 0.31, 0),
S2 = c(0, 0, 0, 4e-04, 0, 0),
S3 = c(0, 0, 0, 0, 0.02, 0),
S4 = c(0.103, 0, 0, 0, 0.038, 0),
S5 = c(0.244, 0.1354, 0.0874, 0.01, 0.003, 0.008),
S6 = c(0.001, 0.094, 0.139, 0.016, 0.03, 0.001)),
class = "data.frame",
row.names = c("OTU1", "OTU2", "OTU3", "OTU4", "OTU5", "OTU6"))

> environment <- structure(list(category = structure(c(1L, 1L, 3L, 3L, 2L, 2L),
.Label = c("Leaf", "Root", "Soil"), class = "factor"),
season = structure(c(1L, 1L, 1L, 2L, 2L, 1L),
.Label = c("dry", "wet"), class = "factor")),
class = "data.frame",
row.names = c("S1", "S2", "S3", "S4", "S5", "S6"))

> pheatmap(matrix,
annotation_col = environment,
cluster_rows = TRUE,
cluster_cols = FALSE,
show_colnames = TRUE,
cluster_distance_rows = "euclidean",
clustering_method = "complete",
fontsize = 10,
fontsize_row = 10)

![Heatmap][1]

[1]: /media/images/07165394-b1d7-4707-8252-65bfd03a
Question: I am trying to understand what an objective function really is. Every sequence alignment method should have an objective function. When it comes to the needleman-wunsch algorithm, what is precisely the objective function? I know there is a dynamic-programming matrix that should be filled. At each cell (i,j) in the matrix, the value is the maximum of three things: value in cell (i,j-1) - gap penalty, value in cell (i-1,j) -gap penalty, value in cell (i-1,j-1)+substitution score. Is this what is called the objective function ?

Answer: I would consider the objective function to be the resulting scores of the paths through the dynamic-programming matrix. This is what's actually being optimized, after all. That also nicely matches the definition of an objective function, since each path score represents multiple comparison events summarized in a single real value.
Question: Hiiii....If anyone can please helo me how to convert fasta files into fastq using ubuntu....i shall be very thank ful....

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

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

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

[1]: https://sourceforge.net/projects/bbmap/
Question: I would like to have the biotypes (like in ensembl transcripts) but for the NCBI refSeq downloaded from UCSC mysql. Which table should I join and which column should I use?

My current query without the biotype is:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -u genome -D hg19 -N -A -e 'select bin,name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,name2 from refGene'

Answer: Technically the biotype is for a transcript, and not a gene. While in many cases the biotype of all transcripts for a gene will be the same, you get a few that aren't. That said I'm not sure off the top of my head if you can do a join within UCSC between the two tables. You might need to output two datasets (refseq and ensembl) and cross-correlate them yourself with a little script. You can output alternative IDs in both tables, and use whichever you prefer to link the two tables.

**Update**: This is a relatively step-by-step guide of how I generated a table linking ensembl Transcript IDs to RefSeq IDs:

UCSC Table Browser: https://genome.ucsc.edu/cgi-bin/hgTables

**Group**: Gene and Gene Predictions **Group**: Ensembl **Table**: ensGene

**Output format**: selected fields from primary and related tables

1. Click on get output
2. On next page under linked tables click the box beside ccdsInfo and then click allow selection from checked tables
3. More tables come up for linking click on ccds id under the CCDS table info fields. Also click on the table knownTo refseq and allow selection from checked tables again
4. Under known to refseq click both fields (primary id and value)
5. Click get output near the top of the page under the fields for the ensembl table

You'll end up with ensembl transcript IDs in the first column and a list of NM IDs in the final column. You can then process this however you like to get a mapping of refseq IDs to Ensembl. I didn't poke around enough to see if I could find a linked table to get biotype IDs, that may be easier to get from BioMart on the ensembl website itself. With those two files you should be able to parse them as tab delimited data and create a mapping file, associate biotypes, etc.with a fairly simple perl/python/scripting language of choice script.
Question: There are many experiments where samples are stored in 384-well plates (or 96-well plates etc.) and a measurement is made on each sample. Typical technical biases may affect spatial regions of the plate, for instance the borders (evaporation, ...), or large bands when the thermocycling (in PCR) is preformed by separate independent controllers, etc. In some more complicated experiments, each well may contain a NGS library, for which multiple quality control scores can be computed.

I searched for `R` functions for pretty printing of values arranged in the geometry of a multiwell plate and I was surprised that I did not find anything obvious. Perhaps because the solution is too trivial ? Still, I would be happy to use a well-thought packaged function instead of a quick-and-dirty ad-hoc home-made wheel reinvention.

Here is an example below. The problems with this script are that 1) it requires flipping some levels that should stay as they are otherwise and 2) SVG or PDF version had a very ugly appearance as each square was strongly smoothed...

*My question is whether there are good packaged tools (in R) to do such plots and more operations on representations of multiwell plates.*

plate <- data.frame( Row = rep(LETTERS[1:16], 24)
, Col = unlist(lapply(1:24, rep, 16)))

plate$Row <- factor(plate$Row, levels=rev(levels(plate$Row)))

set.seed(1)
plate[001:096, "val"] <- rpois(96,10)
plate[097:384, "val"] <- rpois(96, 6)

png("plate.png", h=400,w=600)
ggplot2::ggplot(plate, ggplot2::aes(Col, Row, fill=val)) + ggplot2::geom_raster()
dev.off()

![384-well plate with random values][1]

[1]: https://gist.githubusercontent.com/charles-plessy/6325a878429fef8aae6d9ed25f0054ef/raw/f56cb303524f8c943d86a62583cadea17ea4959f/plate.png

Answer: I made an R package for plotting platemaps with a few utility functions. Mainly just a wrapper round some ggplot2 code, with a few functions for converting between formats (matrix <-> Well_ID+value).

https://github.com/swarchal/platetools
Question: I would like to extract VCF subset with only those variants that are present in one particular sample. E.g. i have three sample VCF with ten variants. Only one variant has been detected in the first sample, so i would like to subset VCF so that only this variant is kept.
Pseudo code would be something like:

bcftools query -f'%CHROM:%POS %INFO [%GT]\n' -i'GT="alt in sample1"' file.vcf

The desired output is VCF with three samples but only variants present in one of the samples, in this case `sample1`

Answer: Alternative solution is to use `bcftools +split` plugin, e.g.:

bcftools +split file.vcf -Ob -o testDir -i'GT="alt"'
Question: So, I have a directory:

Pipeline

- database: contains the custom blast database
- scripts
- genomes
- output

I want my local blast.sh script ( in the script folder) to take all the *.fna of the genomes folder and put results in the output folder.

So i specified the following:

for F in ../genomes/*.fna
do
blastn -query $F -db ../database/blastdb/dbGOI.fasta -outfmt 5 -out ../output/${F%.*}.xml
done

But when I run this, all of my output .xml files are in the genome folder and I do not know why since I specified that in the code.

Am I missing something?

Answer: So, since that seemed to be a bug, I used the following workaround:

blastn -query $F -db ../database/blastdb/dbGOI.fasta -outfmt 5 >> ../output/results.xml

since I used cat to put all theoutput files together, I save a line of code.
Question: I want to have a bigwig file converted into gff after having the genome coordinates and genome annotation files converted from Hg19 to Hg18. One way, would be to convert bigwig into BED, use the liftover tool and then have the BED file converted into gff. Could someone please suggest a method to convert bigwig to BED, nothing available on the net.

Answer: bigWigToBedGraph might be what you're looking for.

bedgraph format
Question: I'm writing my RNA-seq pipeline by Snakemake.when I was writing the last part"rule fpkm" which is calculate fpkm value from bam files. I got the error:

MissingInputException in line 3 of /root/s/r/snakemake/my_rnaseq_data/Snakefile:
Missing input files for rule all:
05_ft/wt2_transcript.gtf
05_ft/wt1_transcript.gtf
05_ft/wt2_gene.gtf
05_ft/epcr1_gene.gtf
05_ft/wt1_gene.gtf
05_ft/epcr2_transcript.gtf
05_ft/epcr1_transcript.gtf
05_ft/epcr2_gene.gtf

Here is my Snakefile:

SBT=["wt1","wt2","epcr1","epcr2"]

rule all:
input:
expand("02_clean/{nico}_1.paired.fq", nico=SBT),
expand("02_clean/{nico}_2.paired.fq", nico=SBT),
expand("03_align/{nico}.bam", nico=SBT),
expand("04_exp/{nico}_count.txt", nico=SBT),
expand("05_ft/{nico}_gene.gtf", nico=SBT),
expand("05_ft/{nico}_transcript.gtf", nico=SBT)

rule trim:
input:
"01_raw/{nico}_1.fastq",
"01_raw/{nico}_2.fastq"
output:
"02_clean/{nico}_1.paired.fq.gz",
"02_clean/{nico}_1.unpaired.fq.gz",
"02_clean/{nico}_2.paired.fq.gz",
"02_clean/{nico}_2.unpaired.fq.gz",
shell:
"java -jar /software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 16 {input[0]} {input[1]} {output[0]} {output[1]} {output[2]} {output[3]} ILLUMINACLIP:/software/Trimmomatic-0.36/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 &"

rule gzip:
input:
"02_clean/{nico}_1.paired.fq.gz",
"02_clean/{nico}_2.paired.fq.gz"
output:
"02_clean/{nico}_1.paired.fq",
"02_clean/{nico}_2.paired.fq"
run:
shell("gzip -d {input[0]} > {output[0]}")
shell("gzip -d {input[1]} > {output[1]}")

rule map:
input:
"02_clean/{nico}_1.paired.fq",
"02_clean/{nico}_2.paired.fq"
output:
"03_align/{nico}.sam"
log:
"logs/map/{nico}.log"
threads: 40
shell:
"hisat2 -p 20 --dta -x /root/s/r/p/A_th/WT-Al_VS_WT-CK/index/tair10 -1 {input[0]} -2 {input[1]} -S {output} >{log} 2>&1 &"

rule sort2bam:
input:
"03_align/{nico}.sam"
output:
"03_align/{nico}.bam"
threads:30
shell:
"samtools sort -@ 20 -m 20G -o {output} {input} "

rule count:
input:
"03_align/{nico}.bam"
output:
"04_exp/{nico}_count.txt"
shell:
"featureCounts -T 10 -p -t exon -g gene_id -a /root/s/r/p/A_th/WT-Al_VS_WT-CK/genome/tair10.gtf -o {output} {input}"

rule fpkm:
input:
"03_align/{nico}.bam"
output:
"05_ft/{nico}_gene.gtf"
"05_ft/{nico}_transcript.gtf"
shell:
"stringtie -e -p 30 -G /root/s/r/p/A_th/WT-Al_VS_WT-CK/index/tair10 -A {output[0]} -o {output[1]} {input}"

Here is my partial directory:

|-- 03_align
| |-- epcr1.bam
| |-- epcr1.sam
| |-- epcr2.bam
| |-- epcr2.sam
| |-- wt1.bam
| |-- wt1.sam
| |-- wt2.bam
| `-- wt2.sam
|-- 04_exp

And the bam files are existed as I run the Snakefile before I add the 'rule fpkm' part.
So before the last part ,it goes fine. But when I add 'rule fpkm' it report error like I describe above.Any ideals are welcomed. Thank you!

Answer: I know the solution, add `,` at the tail of `05_ft/{nico}_gene.gtf`
Question: Is there a ready-made command-line program that can output sequence lengths from a fasta file? Like so:

```
YKR054C,4092
YLR106C,4910
...
```

I've seen [examples][1] of people trying to roll their own using `awk`, but I would rather not do that if I can avoid it. Between Emboss, Biopython, Bioperl, etc there must be something that can do this, right? What I'm looking for would be the equivalent of what `blast_formatter` can do for blast hits, something that can extract simple information that's already there in a format specified by user. Thank you.

[1]: http://stackoverflow.com/questions/23992646/sequence-length-of-fasta-file

Answer: Heng Li's [bioawk][1] might be a good compromise between the simplicity of awk and readability of programming languages.

bioawk -c fastx '{ print $name, length($seq) }' < sequences.fa

[1]: https://github.com/lh3/bioawk
Question: I am supposed to analyze some metabarcoding reads. However, the forward and reverse reads are unable to be merged due to lack of overlap. I was informed that this is because the sequence was too long, so the forward and reverse sequences couldn't extend far enough to overlap adequately.
My question is, what would be the problem of using the forward reads alone, as if they were single-end? My first thought is that they are too short to be a legitimate barcode sequence for identifying taxa. But I'm not sure. It's COI. I gather from Meusnier et al. 2008 that a 95% success rate of species identification was obtained with 250-bp mini barcodes. My forward sequences are that long. But this region by itself has not been tested for specificity. How does this impact its reliability?
Thank you for your input.

Answer: Using just the forward read is a good idea. Just watch out that depending on your library preparation method read 1 might not correspond to the forward direction but forward and reverse is mixes ~50:50. Here trimming the forward primer on read one and two, using e.g. Cutadapt and then reverse complementing the read 2 can help.
If you are concerned about the reliability of the identification with just one direction, you could also analyze read 2 the same way and compare results, or fill in the missing basses in between the reads with Ns to obtain a "full-length sequence". If you do so, however, make sure to apply strict filtering afterward to discard reads of poor quality especially read 2 ends. You can also concatenate sequences, and "reformat" the sequences in your reference database to match these. But this is maybe a bit much effort, using only forward direction should be sufficient in most cases my opinion =)
Question: I have five exomes(in bam format), each a different sample, that I would like to merge into one bam. None of the bam have @RG's in the headers, I am assuming I must add this information first? I also assume this is how I can pull apart the variant information after I processes it with GATK? The Picard documentation states you must pass in a value for UNMAPPED_BAM (definition Original SAM or BAM file of unmapped reads, which must be in queryname order. Required.) I don't undstand what this file is? When I did my alignment should I have dumped the unaligned reads to a file? If so can merge them in a way that does not require unaligned reads to be present.

Answer: I've successfully used the following to merge multiple bam files:

samtools merge merged.bam in.1.bam in.2.bam ...

Although they had @SQs, none had @RGs. Do you have no headers or just no @RGs?
Question: Hello,

I am reading the STAR manual and I came across the publication for STAR-fusion. I could not understand how different they are from each other. Could someone help me understand the difference?

Thank you in advance.

Answer: STAR aligner is a splice aware aligner that aligns junction reads using a guide GTF. [STAR-fusion][1] allows identification of [fusion transcripts][2]. In many cases you get fusion transcripts due to chromosomal translocation (cancer) or trans-splicing etc. Other splice aware aligners ([Tophat-Fusion][3]) also have fusion engine for same purpose.

[1]: https://github.com/STAR-Fusion/STAR-Fusion/wiki
[2]: https://en.wikipedia.org/wiki/Fusion_transcript
[3]: http://ccb.jhu.edu/software/tophat/fusion_index.shtml
Question: Hi,

Hope someone may be able to point me in the right direction!

I've generated some polygenic risk scores for schizophrenia in a cohort of individuals with borderline intellectual disability. Trying to figure out a way to validate the scores I'm getting out of PRSice. Is there a way to do this?

My current results (with --all-score) look something like this:
https://ibb.co/jv9YQP5
(sorry couldn't place the table here in an easily legible format)

There are some patterns in my results, with some individuals having identical scores at particular p-values (is this just due to the stats involved in generating the scores?), and the scores are pretty low (clustering around 0) with some changing from positive to negative depending on p-value threshold. Is this typical of polygenic risk scores?

I've followed Marees et al 2018 (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6001694/) and Sam Choi's tutorail (https://choishingwan.github.io/PRS-Tutorial/) to generate this, but I'm thinking I may have made errors somewhere along the line and would like to check my results somehow!

Any advice would be much appreciated!

Answer: Your sample size is simply too small. That will usually not give you any reasonable result. I suspect that as you've less than 100 samples, the clumping LD varies quite a lot, lead to more variants being removed than normally would (which lead to the small variant count). One way to get around that will be to use a reference panel (e.g. 1000 genome) which has a larger sample size. Though at the end of the day, with only 100 samples, I don't expect you to get any meaningful results.

Given there are only 265 variants after clumping, it is very possible that you'd only have one or two SNP included in a particular threshold, which can then lead to the identical PRS. So I am not too surprise about your results.
Question: we did a 16S RNA analysis using MiSeq kit version 3. This kit makes paired-end reads with a size of 300pb. By merging both reads, we can cover 535 pb of region V3-V5.

As you can see in the following plot, qualities of reads are really bad after the nucleotid 250.
Is it normal ? Perhaps the strategy prefers reads lenghs to the quality ? Did you get something similiar with the same kit ?

![V3-V5 region][2]
![fastqt quality plot][1]
After merging with Flash
![After merging][3]

[1]: https://framapic.org/7luICTwCp8FU/dyD2VMKC7cXW.png
[2]: https://framapic.org/TASsvk65DaIc/xpkU8WlSSVWX.png
[3]: https://framapic.org/JNjXjmVbVhHm/1oVAPSZiYRSu.png

Answer: As pointed in the comments, Illumina's problems with MiSeq v3 kit and 2x300bp are known and old, see for example [here][1] and [here][2]. Informally, (at least some) Illumina technicians will tell you to stay with 2x250bp, I am not sure if someone at Illumina said so on record.

[1]: http://blog.mothur.org/2014/09/11/Why-such-a-large-distance-matrix/
[2]: http://omicsomics.blogspot.com.br/2015/12/loose-2015-threads-1-miseq-2x300-issues.html