cannin/biostars_qa huggingface dataset:
Question: Hi,

I was hoping someone with nextflow experience could help me with this issue.

My script is a run of the mill Hisat2/Stringtie nextflow script. However my issue is that the final process of the script takes one file at random from the previous channel, and finishes the script without any errors. Here is a portion of the code showing where the inputs to the final process originate from:

#!/usr/bin/env nextflow

params.genome = "Reference/chr22.fa"
genome_fasta = files( params.genome )

params.annot = "Annotation/chr22.gtf"
Channel
.fromPath( params.annot )
.into { gtf1; gtf2; gtf3 }

params.reads = "trimmed_reads/*_r{1,2}.trimmed.fastq.gz"
Channel
.fromFilePairs( params.reads )
.set { read_ch }

......
extract exons, splice sites, align with hisat2, pipe to bam
......

process Sort_Index_Bams {
publishDir "BAMS/", mode:'copy'

input:
set val(key), file(bam) from hisat_bams

output:
set val(key), file("${key}.bam") into hisat_bams1
file "${key}.bam.bai" into indexed

script:
def avail_mem = task.memory == null ? '' : "-m ${task.memory.toBytes() / task.cpus}"
"""
samtools sort \\
$bam \\
-@ ${task.cpus} $avail_mem \\
-o ${key}.bam
samtools index ${key}.bam
"""
}

hisat_bams1.into { hisat_bams2; hisat_bams3 }

process Assemble_Transcripts{
publishDir "Assembly/", mode:'copy'

input:
set val(key), file(bam) from hisat_bams2
file(gtf) from gtf2

output:
file("${key}.gtf") into hisat_transcripts

script:
"""
stringtie \
${bam} \
-G ${gtf} \
-l ${key} \
-o ${key}.gtf \
-p ${task.cpus} \
"""
}

I have tried to alter the final process as such:

- ```set val(key), file (bam) from hisat_bams2.collect()```. This returned an error of "Input tuple does not match input set cardinality declared by process Assemble_Transcripts".
- Set the input GTF file to a value channel as described on this stackoverflow post [here][1]. The user reported same issue, however it did not solve my problem.

Any suggestions would be greatly, greatly appreciated.

Regards,

Barry

[1]: https://stackoverflow.com/questions/57535349/issue-with-bwa-mem-process-not-running-on-all-output-files-from-previous-process

Answer: instead of

set val(key), file(bam) from hisat_bams2
file(gtf) from gtf2

I would go for

set key,bam,gtf from hisat_bams2.combine(gtf2)

because here, gtf2 is a Channel, not a file.
Question: Hi there,

In variant annotating with SnpEff + dbNSFP, there are many variants assigned with more than one value in some fields. For example, as highlighted in yellow as below, SIFT_Pred has 5 "T" (but only one SIFT_Score), Polyphen2_HDIV_pred has 2 "B".

If I understand it correctly, each variant should only have one SIFT_score, SIFT_Pred, etc. Did I miss anything?

Thanks in advance.

7 140301731 . G T 2918.77 PASS AC=2;AF=1.00;AN=2;DP=95;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQ0=0;POSITIVE_TRAIN_SITE;QD=30.72;VQSLOD=6.96;culprit=FS;ANN=T|missense_variant|MODERATE|DENND2A|ENSG00000146966|transcript|ENST00000275884|protein_coding|2/19|c.467C>A|p.Pro156His|885/3735|467/3030|156/1009||;dbNSFP_rs_dbSNP141=rs269243;dbNSFP_CADD_phred=11.67;dbNSFP_ESP6500_AA_AF=0.989765;dbNSFP_ExAC_AF=9.549e-01;dbNSFP_UniSNP_ids=rs269243,rs1681775,rs52835101,rs58245328,rs269243;dbNSFP_SIFT_pred=T,T,T,T,T;dbNSFP_Polyphen2_HDIV_rankscore=0.02656;dbNSFP_Polyphen2_HDIV_pred=B,B;dbNSFP_ESP6500_EA_AF=0.936014;dbNSFP_SIFT_converted_rankscore=0.12818;dbNSFP_Polyphen2_HVAR_rankscore=0.01281;dbNSFP_1000Gp1_AF=0.9647435897435898;dbNSFP_Polyphen2_HVAR_pred=B,B GT:AD:DP:GQ:PL 1/1:0,95:95:99:2947,286,0

Answer: dbNSFP provides multiple values for those scores and SnpEff simply annotates them all.

The reason dbNSFP provides different values is related to how each score is calculated. For instance the number of PolyPhen entries usually the same as the number of Uniprot IDs, similarly the number of SIFT entries relates to the number of "Ensembl_transcriptID" entries.
Question: Dear All

I have a problem that we can not solved for a while in installing bioawk. I follow the general pipeline to install bioawk on ubuntu server:

- clone from git hub to the Download directory
- make

After that it is where we are stuck:

When I run the command:

mv awk book

I got error awk does not exist! When I install gawk and try `mv gawk` from `/usr/bin` to `home/arsalan/download/bioawk`, it works but when I `cp bioawk` to `usr/local/bin` and try to run it if fails and since gawk has already moved from `/usr/bin`. By typing `gawk` I got the error it is not installed! and typing `bioawk` I got error command not found!

Many thanks in advance

Answer: Clone bioawk repo, run make and make bioawk executable

```
git clone https://github.com/lh3/bioawk.git
cd bioawk
sudo make
sudo chmod +x bioawk
./bioawk
```
Question: I am trying to loop paired end reads in hisat2 with an output of sorted bam using a pipe.

I am a novice at coding but am keen to improve my skills by learning how to loop.

I have created this loop using the answers from https://www.biostars.org/p/223404/

for f in `ls *.fastq.gz | sed 's/_[12].fastq.gz//g' | sort -u`
do
echo hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 ${f}_1.fastq.gz -2 ${f}_2.fastq.gz | samtools sort -o ${f}.hisat.sorted.bam
done

For some reason the loop does work without the pipe but not with. How can i configure this loop to work with the piping as I do not want to have the intermediate files?

Many thanks in advance

Answer: Remove the `echo` as you want to execute rather than print the command.
If you want "status updates" do

echo $f && hisat2 (...)
Question: I often see the word depth in the manuals of the tools for NGS, what is its meaning ?

thanks.

Answer: It's also worth noting that you can use read depth to infer copy number, given a sufficient number of reads. On average, you'll get 1.5x more reads from a triploid region than you will from the rest of the diploid genome.
Question: Initially I [asked a question][1] of how to make a dotplot from a list of pathways.
I was able to generate the graph, but now I can not make the groups side by side using the pathways of reference.

I do not know if that's how I tabulated my results, or if it's a problem with my code, could anyone help me?

**My data**

**My code**

data <- read.csv("S3.csv", sep =";", header = TRUE, stringsAsFactors = FALSE)
library("ggplot2")
S3<- ggplot(data, aes(x="Groups", y=Pathways, size=DEGs, color=FDR)) + geom_point(alpha = 0.8) + facet_grid(Group ~ .) + theme_classic()
S3
S3 = S3+scale_color_gradient(low = "red2", high = "mediumblue", space = "Lab", limit = c(0.000000000000000000000000000000000000000000000004, 0.03))
S3+scale_size(range = c(2, 8))
S3

**My result**

**How do I want the dotplot look like (Edited image)**

Best regards,
Leite

[1]: https://www.biostars.org/p/359307/#360478

Answer: You need to just use facet_wrap(Group ~ .) instead of facet_grid
Question: Is it possible to skip the estimateSizeFactor step in Deseq2 and define your own size factors? If I just replace the size factor vector with my own values, will the down-stream DE steps use these values? I am pretty new with R so it would be very helpful if I can see some example code on how to do this.

Answer: Yes you can, that is what I can read from the below link

Take a look [here][1] and [here][2]

There was a [post][3] from where Damian you already replied few years back where I can see Jeremy replied that you can set the size factor manually as well with `sizeFactors`

*dds <- estimatesizefactors(dds)<="" code=""> See DESeq for a description of the use of size factors in the GLM. One should call this function after DESeqDataSet unless size factors are manually specified with sizeFactors. Alternatively, gene-specific normalization factors for each sample can be provided using normalizationFactors which will always preempt sizeFactors in calculations. Internally, the function calls estimateSizeFactorsForMatrix, which provides more details on the calculation.*

[1]: http://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/estimateSizeFactors
[2]: http://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/sizeFactors
[3]: https://www.biostars.org/p/12273/
Question: I am using HOMER to get read mapping around TE (well, I need to make a heatmap). I am using following command:

annotatePeaks.pl XXX.bed XXX(genome) -bedGraph XXX.bedgraph -size 200 -hist 2 -ghist -noadj -fragLength 0 > XXX.txt

Before running this command, I converted my bed file using changeNewLine.pl and configured my genome in HOMER. I have been using this tool and so this command for couple of weeks now. From terminal result it seems ok. Here is the result in terminal:

Peak file = XXX.bed
Genome = XXX
Organism = null
bedGraph Files:
XXX.bedgraph
Peak Region set to 200
-----------------------------------------------------
Histogram mode activated (bin size = 2 bp)
-----------------------------------------------------
Will create histogram for each gene
Will NOT normalize tag counts
Fragment Length set to 0
Peak/BED file conversion summary:
BED/Header formatted lines: 4203
peakfile formatted lines: 0
Duplicated Peak IDs: 0

Peak File Statistics:
Total Peaks: 4203
Redundant Peak IDs: 0
Peaks lacking information: 0 (need at least 5 columns per peak)
Peaks with misformatted coordinates: 0 (should be integer)
Peaks with misformatted strand: 0 (should be either +/- or 0/1)

Peak file looks good!

Resizing peaks...
Reading Positions...
-----------------------
Compiling per bp Histograms...
Finding Tags in Peaks from each directory...

Order of experiments in output file:
XXX.bedgraph

But, the number of rows in my output file and number in input bed file are not same. Number in output file is much less. Would you please tell me what went wrong and how can I get all rows in my output file?

Thank you very much and I am very sorry to bother.

Answer: If you just want to get gene annotation for your BED regions, you only need to use annotatePeaks like this:

perl annotatePeaks.pl MyBED.bed mm9 -annStats AnnotationStats.txt > Annotation.tsv

This always produces an equal number of annotated regions as per my BED file. The correct genome has to be used of course.

HOMER will assume that the BED file regions define your peaks. I don't believe that you have t supply all of the other information such as the bedGraph, peak size, etc.
Question: Hello, I have a bed file and need to merge with strand specificity, but I am not getting an output. My command is:

bedtools merge -i input.bed -s > merged.bed

My bed file looks like :

```
chr1 65702 65714 +
chr1 71879 71891 -
chr1 97442 97454 -
chr1 97552 97564 -
chr1 125938 125950 +
chr1 134489 134501 -
chr1 165663 165675 -
chr1 173419 173431 +
chr1 176257 176269 +
chr1 176263 176275 +
```

Answer: Try this:

awk '{OFS="\t"; print $1,$2,$3,".",".",$4}' INPUT.bed | bedtools merge -i - -s

[Bed format][1] contains strand information in column 6.

[1]: http://genome.ucsc.edu/FAQ/FAQformat.html#format1
Question: Dear all,

I'm analyzing eukaryotic metatranscriptomics HiSeq (2x100 bp) data from soil (the cDNAs are from eukaryotic poly-A mRNAs). The FastQC report of the R1 reads shows that there is a sudden sequence quality drop around the 50th nucleotide, even some unknown bases (Ns), and high Kmer content (see [here][1]). However, this drop in quality and unknown bases are absent in the R2 reads ([there][2]).

I wonder if it is normal to have something like this, and if it is not, what should I do to improve the read quality?

[1]: https://drive.google.com/file/d/0B1HJ6hkNttmvY3VLR1A4SlBVSk0/view?usp=sharing
[2]: https://drive.google.com/file/d/0B1HJ6hkNttmvU05HbTNrczFOSXc/view?usp=sharing

Answer: This is likely due to a bubble going through as the image was being acquired. This actually used to be fairly common. Whether this will be an issue or not will depend a bit on exactly what you want to do with the data.
Question: I am from a Biology background. To learn command line, GitHub, and Linux I have installed Ubuntu on Windows 10 laptop.

When I run the app message pop-up

To run a command as administrator (user "root"), use "sudo ".
See "man sudo_root" for details.

name@LAPTOP-xxxxx:

when I look the location running `~$ pwd`

it shows `/home/name`

1. Using the ubuntu app — how can I move to D drive and different folders in D drive?

2. Also, what is the difference between running Ubuntu on Windows and **Windows Command Prompt** (cmd.exe)? Are they similar?

3. Can I run command line script on Ubuntu and cmd?

I tried to Google these. however unable to find a friendly answer.

I would be grateful if someone shares what is happening when running Ubuntu on Windows and any relevant tutorial if possible.

Answer: Windows drives are mounted under `/mnt`, for example:

ls /mnt/c

The Linux subsystem can be thought of as a separate, simulated Linux environment.

You cannot run Linux commands on Windows. But because the Linux is simulated within Windows, you can start some Windows programs from it. For example typing:

explorer.exe .

in the Linux bash, will open a Windows Explorer instance that shows the current **linux** directory. Very handy shortcut.

Or, if you have Visual Studio Code you can start it from the Linux Subsystem just by typing `code`. The newly opened Visual Code editor runs in Windows, but it is connected to Linux, code that is executed from the editor will run under Linux.

In general you should think of the two as separate environments that only share files.
Question: Hello. I've been recently trying to obtain experimental data on protein interactions, up to the pdb entry. I've come to realize, there might not be such an easy way of obtaining this data - For example Stringdb or intAct only provide Uniprot or some other accession entries (If I'm not mistaken), And obtaining exact pdb from just the Uniprot is currently not possible - at least not precisely.

Is there any other way of obtaining experimental PPI pdb data?

Thank you very much in advance.

Answer: I am not sure what you mean by not possible to get exact PDB from UniProt. PDB is searchable by Uniprot entry. UniProt has a PDB cross-reference file and PDB also provides a list of structures to Uniprot mapping.

Also Andrew Martin's group at UCL provides a web interface to the PDB/UniProt mapping with REST access.

You may also be interested in the EBI's SIFT (Structure integration with function, taxonomy and sequence) resource.
Question: Hello guys,

It may seem like a basic question, but this is causing confusion.

What is the difference between fold change and Log fold change?

Regards,

Answer: To follow up on the value added by the symmetry in log2-transformed fold changes, you might also look at a [volcano plot][1] of features, a common way to represent p-values vs log2-fold change.

For example:

![enter image description here][2]

Symmetry makes visual interpretation of this type of plot — comparison and labeling of features with significant fold change
and p-value thresholds — much simpler. One can see immediately which features are both statistically significant and likely to be significantly up- or down-regulated.

[1]: https://en.wikipedia.org/wiki/Volcano_plot_(statistics)
[2]: https://image.ibb.co/fC2B4S/volcano.png
Question: Dear all,

I'm struggling in finding how to obtain information about COSMIC database after I annotate a vcf using VEP.
The command that I used is:

```
vep	-i Mutect2_unfiltered_10643_vs_2434.vcf.gz
-o Mutect2_unfiltered_10643_vs_2434_VEP.ann.vcf
--assembly GRCh38
--species homo_sapiens
--offline --cache --cache_version 99
--dir_cache /.vep --everything --filter_common
--fork 4 --format vcf --per_gene --stats_file Mutect2_unfiltered_10643_vs_2434_VEP.summary.html
--total_length --vcf
```

I get back a vcf annotated file, but I'm not able to find any COSMIC information. Maybe there is a way to filter only for COSMIC?

Thank you in advance,
Youssef

Answer: The COSMIC variants should appear listed as SOMATIC under the co-located variants. You may find that your `--filter_common` option is excluding some COSMIC variants. I would recommend running without this then using the [VEP filter tool](http://www.ensembl.org/info/docs/tools/vep/script/vep_filter.html) with `--filter "SOMATIC"` to find COSMIC variants.
Question: I have a dataframe in R that looks like this:

V1 T1 T2 T3 T4 T5
CXCL6 0.8536601 1.0903336 3.7633042 5.5800459 5.8477150
PPBP 0.7739450 0.3587961 0.5073359 0.2743522 0.6221722
CXCL10 0.1258370 -0.3535165 -0.7460387 3.5604672 0.1971432
CXCL11 -0.2563139 0.7117200 0.0000000 -0.2288303 0.9955557
CXCL12 0.6181279 1.7529310 1.7637760 1.2752787 1.2284810

I want to keep the rows that have values only between -1 and 1.

I have tried this command but unfortunately it does not work.

condition1 <- Genes[,c(2:6)] <=-1 & Genes[,c(2:6)] >=1
Genes <- Genes[condition1,]

Can someone tell me where I am wrong so that I can successfully filter my dataframe.

Answer: Your condition is not clear, but if you want rows that contain values that are ONLY between -1 and 1, this should work :

Genes[which(rowSums( (Genes[,2:6]>=-1) & (Genes[,2:6]<=1))==5),]
Question: Hello biostars,

I am trying to understand what rho means and how it is calculated in coalescent simulations.
We know for Hundon's 'ms' program for example, that cross-over rate parameter is measured by 4N0 * r (r is probability of cross-over per generation between the ends of the locus).
Now I want to simulate a locus in ms using -r parameter and I know average recombination rate for my species. How can I calculate per base per generation crossover rate?

Average recombination rate is in cM\Mb, so if, let's say, human recombination rate is 2 cM\Mb, genome size is 3 * 10^9 , then 2 / 1 * 10 ^6 = 2 * 10^-6 recombinations per base pair. This is way higher than a number used by default in different papers (such as 2*10^-8). And if I divide it by 3 * 10^9, the number is so small that it does not make sense at all. I know my calculations are incorrect, but I just don't see another way to derive this parameter.

Answer: Hey,

I don't know if you've already figured this out, but aside from the complexities of estimating recombination rate - the big problem in your estimate is conflating cM/Mb with crossovers/Mb. The genetic distance of 1 centimorgan represents only a 1% chance of crossover in the region. If you take 2cM/Mb to be (2 x 10^-2 crossovers) / (1 x 10^-6 bp), you get the 2 x 10^-8 crossovers/bp.
Question: Hi,

Any idea if there's a way to get the Instrument type (i.e. NextSeq, HiSeq, etc), from the instrument name field in the Fastq read header? - I've tried looking for a map of sorts in Illumina's documentation, but to no avail.

Thanks,

Answer: There was a post @SeqAnswers where we had a recap of the types of sequencer specific fastq headers (my google-fu has not turned up that post yet). Following is a very rough approximation of start of fastq headers that are sequencer specific.

Edit: Just want to add that it is possible to change the default read header in one of the set up files (so the following may not always be true).

@HWI-Mxxxx or @Mxxxx - MiSeq

@HWUSI - GAIIx

@HWI-Dxxxx - HiSeq 2000/2500

@Kxxxx - HiSeq 3000(?)/4000

@Nxxxx - NextSeq 500/550

@Axxxxx - NovaSeq

@Vxxxxx = NextSeq 2000

Edit (08/2022)

@AAxxxxx - NextSeq 2000 P1/P2/P3

@Hxxxxxx - NovaSeq S1/S2/S4
Question: **Intro:**
Yeast chromosome 12 contains a really long region of rDNA. There are about 150-200 or greater copies of a 9.1kb repeat, so about 1.5M bases in total, which is the majority of the chromosome. The published yeast reference genome has only 2 copies of the repeat.

**My situation:**
I just did whole genome sequencing of some yeast clones. I have on about 70X coverage in most of the genome. If you figure that there are 200 copies of the rDNA repeat on chromosome 12, but only 2 copies represented in the reference sequence, then in this region I should see a read depth of roughly 70*100 = 7000. And in fact, this is pretty much what I see, here's a plot of read depth in that area. To generate this plot I mapped reads with `bwa` and then used `samtools depth` and then plotted in R.

![enter image description here][1]

**My question:**
The plot looks obviously artificially cut off at a read depth of ~8000--why? I didn't specify any limits in creating the plot. **So where is the cutoff happening?** Is it something samtools depth is doing? Something bwa did? Or something about how the reads are generated during Illumina sequencing?

Last thing to mentions: the highest value is not exactly 8000, there are a few positions with depths of 8002, just to make things weirder. And I've shown just one example, but it happens on to all clones I sequenced.

Thanks for any help!

[1]: https://i.imgur.com/I2ReY5E.jpg

Answer: Usage: samtools depth [options] in1.bam [in2.bam [...]]
Options:
(...)

-d/-m maximum coverage depth [8000]. If 0, depth is set to the maximum
integer value, effectively removing any depth limit.
(...)
Question: Given an indexed BAM file, is the following command guaranteed to give the full and correct listing of (UCSC) chromosome names for reads in the BAM file?

$ samtools view -H foo.bam | cut -f2 | grep '^SN:chr' | sed s'/SN://'

I'm trying to think of any gotchas or custom implementations. Are there other tags that could be used for chromosome names in BAM file headers, or other issues I'm not thinking of?

Answer: If the file is already indexed, you can use the `idxstats` tool.

samtools idxstats NA18152.bam | head -n 3
chr1 247249719 123894 0
chr2 242951149 97215 0
chr3 199501827 81334 0

Then:

samtools idxstats NA18152.bam | cut -f 1 | head -3
chr1
chr2
chr3
Question: Hi everyone, i am using bioperl to do blast 2 sequence. I need to output as a tabular format. I saw some post where they mention we need to use -m8 to get tabular output. But in my code its not working. Can anyone help me out what i am missing in this code.

```
use Bio::Tools::Run::StandAloneBlastPlus;
#use Excel::Writer::XLSX;
use strict;
use Bio::SearchIO;

my $factory = Bio::Tools::Run::StandAloneBlastPlus->new(
-db_data =>'test.fasta', # Specifies the file to use for database
-create =>1 # Creates a new temporary database
);
my $result = $factory->blastn(
-query =>'H1.fasta', # Specifies the query file
-outfmt => -m8,
-outfile =>'out_H1' # Specifies the output location
);
$factory->cleanup; # Deletes the temporary database files

```

Answer: As far as I remember the `-outfmt 6` stands for tabular output. See https://www.biostars.org/p/88944/
Question: I am rookie to using limma package. My data looks like this in excel: first column is the gene name, first raw is sample name, and under each sample name is the gene expression value (4 replicates per sample, I have sample wt and mu). I want to find out the significant expressed genes between mu vs wt.

and I write code like this:

sample=read.csv("sample.csv",header=T,row.names=1)
logsample=log2(sample)
design=model.matrix(~0+c(rep('wt',4),rep('mu',4)))
colnames(design)=c("wt","mu")
cm=makeContrasts(mu-wt,levels=design)
fit=lmFit(logsample,design)
fit2=contrasts.fit(fit,cm)
fit3=eBayes(fit2)
result=topTable(fit3,number=Inf,adjust="BH",sort.by="none")

I am not sure I did it correctly. Please help me to check, many thanks!! I log2 transformed the intensity. The data is already normalized by other software; I am not sure I should use makeContrasts or not.

Thanks!!

Answer: Your code looks good to me to find genes differentially expressed between the 4 WT and 4 mu, which I guess it is what you want.

As Devon Ryan suggested, you can leave the intercept and skip the contrasts. This should be equivalent to your code:

```
design=model.matrix(~1+c(rep('wt',4),rep('mu',4)))
colnames(design)=c("wt","mu")
fit=lmFit(logsample,design)
fit<- eBayes(fit)
results<- topTable(fit,number=Inf,adjust="BH",sort.by="none", coef= 2)
```

However, I don't see any harm in explicitly fitting contrasts, actually I find it clearer.

I think in general it's good to look at a couple of genes from the output of `topTable()` and make sure the model is doing what you expect by comparing to the values in the normalized matrix. E.g. in your case:

results$logFC[1]

Should be equal to (given rounding errors):

mean(logsample[1, 1:4]) - mean(logsample[1, 5:8])