cannin/biostars_qa huggingface dataset:
Question: Hi,

I am very new to `nextflow` (I used to work with `snakemake` in the past). I am trying to create a dummy workflow for understanding the basics of the _pipeline creation_ and the first step at designing my own.

For this I want to unzip my fastq files and create 2 dummy repports from the read length of the R2 read.
I have **2** scripts for now, `main_qc.nf` *(including the workflow)* and `modules_qc.nf` *(with all processes)* shown below:
```
#!/usr/bin/env nextflow
nextflow.enable.dsl=2

// Include in workflow
include {
GUNZIP_FASTQ;
GET_READ_LENGTH
} from "./modules_qc.nf"

// Initial parameters
datadir="/data"
bashdir="/shs"

params.sample = "s1"
params.fastq = "$datadir/fastqs/${params.sample}_2.fastq.gz"

workflow {
// Create
Channel
.fromFile(params.fastq)
.set { chFastqFile }

Channel
.of(params.sample)
.set { samples }

GUNZIP_FASTQ(chFastqFile)
GET_READ_LENGTH(samples, GUNZIP_FASTQ.out)
}
```

```
#!/usr/bin/env nextflow
nextflow.enable.dsl=2

// Unzipping files
process GUNZIP_FASTQ {

input:
path target

output:
path "${target.simpleName}.fastq"

script:
"""
gunzip -d -c ${target} > ${target.simpleName}.fastq
"""
}

// Export read length to file
process GET_READ_LENGTH {
input:
val sample_id
path fastq

output:
path "${sample_id}.readLength.txt"

script:
"""
bash ./shs/readLength.sh ${fastq} ${sample_id}.readLength.txt
"""
}
```

I want to first run the `GUNZIP` process on all fastqs for all sample and then create one dummy repports per sample.
`GUNZIP` processes have to run twice ad much as the other processes.

How should I proced?

Thank you very much

Answer: I would write it the following way (not tested)

#!/usr/bin/env nextflow
nextflow.enable.dsl=2

// input is a tsv file with sample and path/to/fastq
params.sample_fastq = ""

workflow {
sample_fastq_ch = Channel.fromPath(params.sample_fastq).splitCsv(header: false,sep='\t')
gunzip_ch = GUNZIP_FASTQ(sample_fastq_ch)
readlen_ch = GET_READ_LENGTH(gunzip_ch.out)
zip_ch = ZIPALL(readlen_ch.out.collect())
}

process GUNZIP_FASTQ {
input:
tuple val(sample),val(fq)
output:
path("${sample}.fastq"),emit:out
script:
"""
gunzip-c ${fq} > ${sample}.fastq
"""
}

process GET_READ_LENGTH {
input:
tuple val(sample),val(fq)
output:
path("${sample_id}.readLength.txt"),emit:out
script:
"""
bash /full/path/to/shs/readLength.sh ${fq} "${sample}.readLength.txt"
"""
}

process ZIPALL {
input:
val(L)

output:
path("output.zip")
script:
"""
zip -j output.zip ${L.join(" ")}
"""
}

a side note : whatever is "readLength.sh", you shouldn't use a software that requires you to gunzip a fastq...
Question: Hello,

I do Variant Calling with freebayes. For comparison with another vcf from another source I need to split called mnp to snp.

For example freebayes called this:

chr2	21232803	.	TG	CA,CG

But I need this:

chr2	21232803	.	T	C
chr2	21232804	.	G	A

How can I do this?

fin swimmer

Answer: here is the solution to your problem with bcftools and [vt][1] and I added an extra record to check the program: (For some reason latest version of vt is crashing on decompose function. Hence bcftools was used. Otherwise, vt alone is fine for this solution)

$ bcftools norm -m -both test.vcf.gz | vt decompose_blocksub - -o output.vcf

output ( VCF headers are not shown):

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
2	21232803	.	T	C	.	.	OLD_CLUMPED=2:21232803:TG/CA
2	21232803	.	T	C	.	.	OLD_CLUMPED=2:21232803:TG/CG
2	21232804	.	G	A	.	.	OLD_CLUMPED=2:21232803:TG/CA
20	763837	.	C	T	.	.	OLD_CLUMPED=20:763837:CA/TG
20	763838	.	A	G	.	.	OLD_CLUMPED=20:763837:CA/TG

input ( VCF headers are not shown):

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
2	21232803	.	TG	CA,CG
20	763837	.	CA	TG

[1]: https://genome.sph.umich.edu/wiki/Vt
Question: Hey guys,

I'm having a problem when using a mouse dbsnp vcf file for variant annotations with GATK's VariantAnnotator.
This dbsnp vcf file has chromosomes notated as chr1, chr2 .... but my reference follows a different notation 1, 2, 3 ... This produces an error when running GATK, like the following:

Input files dbSNP.vcf and reference have incompatible contigs: No overlapping contigs found.

ERROR dbSNP.vcf contigs = [chr1, chr10, chr11, chr12, chr13, chr13random, chr14, chr15, chr16, chr17, chr17random, chr18, chr19, chr1random, chr2, chr3, chr3random, chr4, chr4random, chr5, chr5random, chr6, chr7, chr7random, chr8, chr8random, chr9, chr9random, chrM, chrUnrandom, chrX, chrXrandom, chrY, chrYrandom]

ERROR reference contigs = [1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 2, 3, 4, 5, 6, 7, 8, 9, MT, X, Y, NT166325, NT166464, NT166452, NT166480, NT166448, NT166458, NT166443, NT166466, NT166476, NT166479, NT166478, NT166474, NT166471, NT166445, NT166465, NT166457, NT166470, NT166454, NT166472, NT166449, NT166481, NT166337, NT166459, NT166456, NT166473, NT166461, NT166475, NT166462, NT166444, NT166453, NT166446, NT166469, NT072868, NT166335, NT166467, NT166283, NT166338, NT166340, NT166442, NT166334, NT166286, NT166451, NT166336, NT166339, NT166290, NT053651, NT166450, NT166447, NT166468, NT166460, NT166477, NT166455, NT166291, NT166463, NT166433, NT166402, NT166327, NT166308, NT166309, NT109319, NT166282, NT166314, NT166303, NT112000, NT110857, NT166280, NT166375, NT166311, NT166307, NT166310, NT166323, NT166437, NT166374, NT166364, NT166439, NT166328, NT166438, NT166389, NT162750, NT166436, NT166372, NT166440, NT166326, NT166342, NT166333, NT166435, NT166434, NT166341, NT166376, NT166387, NT166281, NT166313, NT166380, NT166360, NT166441, NT166359, NT166386, NT166356, NT166357, NT166423, NT166384, NT161879, NT161928, NT166388, NT161919, NT166381, NT166367, NT166392, NT166406, NT166365, NT166379, NT166358, NT161913, NT166378, NT166382, NT161926, NT166345, NT166385, NT165789, NT166368, NT166405, NT166390, NT166373, NT166361, NT166348, NT166369, NT161898, NT166417, NT166410, NT166383, NT166362, NT165754, NT166366, NT166363, NT161868, NT166407, NT165793, NT166352, NT161925, NT166412, NT165792, NT161924, NT166422, NT165795, NT166354, NT166350, NT165796, NT161904, NT166370, NT165798, NT165791, NT161885, NT166424, NT166346, NT165794, NT166377, NT166418, NT161877, NT166351, NT166408, NT166349, NT161906, NT166391, NT161892, NT166415, NT165790, NT166420, NT166353, NT166344, NT166371, NT161895, NT166404, NT166413, NT166419, NT161916, NT166347, NT161875, NT161911, NT161897, NT161866, NT166409, NT161872, NT166403, NT161902, NT166414, NT166416, NT166421, NT161923, NT_161937]

Is there any simple way/script that would allow me to change the chromosome notation on the dbsnp file to the same provided on the reference?

Thanks

Answer: Probably:

perl -pe 's/^([^#])/chr\1/'
Question: Hello,
I have some reads aligned to genome and I want to draw a figure showing how many reads mapped to existing annotated lncRNA/mRNA/intergenic regions, something like the figure bellow:

![enter image description here][1]

[1]: /media/images/2a559bae-13ab-4913-aba0-a4c118d5

but a little different that I want the conjunctions to be two-way, like some circos plot but between two set of vertical segments.
Any one have ideas? Thanks!

Answer: Sankey Diagram https://www.r-graph-gallery.com/sankey-diagram.html
Question: Hi,

I aligned my RNA-Seq reads using STAR on human hg19 genome index build with ENSEMBL annotation (`--sjdbGTFfile` option at the genome index generation). So in my SJ.out.tab file I will have information about the number of reads supporting splice junction as described in STAR manual :

```
column 1: chromosome
column 2: first base of the intron (1-based)
column 3: last base of the intron (1-based)
column 4: strand (0: undefined, 1: +, 2: -)
column 5: intron motif: 0: non-canonical; 1: GT/AG, 2: CT/AC, 3: GC/AG, 4: CT/GC, 5:AT/AC, 6: GT/AT
column 6: 0: unannotated, 1: annotated (only if splice junctions database is used)
column 7: number of uniquely mapping reads crossing the junction
column 8: number of multi-mapping reads crossing the junction
column 9: maximum spliced alignment overhang
```

In ENSEMBL annotation, there are all types of RNA (mRNA, tRNA, lncRNA,...). But for my research I've to compute the number of uniquely mapped reads crossing only protein-coding genes splice junction. Is there a way to download splice acceptor/donor position in hg19 for protein-coding genes in order to subtract the read count from the `SJ.out.tab` file?

Thanks

Answer: HISAT actually comes with a python script that can extract splice sites from a GTF file. So just take your GTF file, extract only the types of entries that you want to include, run `extract_splice_sites.py` on the results and then intersect the results appropriately.
Question: I need a way to assign a number to genes representing how broken it is depending on how many SNPs with minor alleles appear on the gene. Based on just a handful of SNPs though, not every SNP.

So lets say for example I know a few SNPs for COMT:

```
Rs165599 AA
Rs6269 AA
Rs737865 AA
Rs4633 TT
Rs769224 GG
Rs4680 AA
```

and some SNPs for STAT3:

```
Rs744166 GG
Rs8069645 AG
Rs9891119 CC
Rs6503691 CC
```

and I wanna find out which gene is more defective (as in which has the most diminshed ability to encode the corresponding protein), how can I do that? These genes obviously have way more SNPs than just those ones, and in every gene I will have many minor alleles so I'm guessing its only certain SNPs that cause problems for the gene. So could I find a number that tells me how much a bad allele of a particular SNP breaks the gene and then use that number to get an overall idea of how malfunctioning the gene is?

Answer: 1. You're looking for [VEP][1] and/or [Clinvar][2]
2) Be careful with your terminology. Evolutionary Biologists might argue SNPs do not break genes - they create them. But more importantly, people who carry variants (with or without a deleterious phenotype) don't like being called "broken", "bad", "defective", "diminished", "problematic" or "malfunctioning" ;-) But before you blush, don't sweat it - we're taught to think this way in school and with good intentions, but it lands you in awkward situations like a presentation I once sat through re. how "defective individuals" carry the Sickle Cell gene, and a girl in the class had to leave because her mother suffers (she was therefor a carrier) and the whole thing was, rightly, upsetting.

[1]: http://www.ensembl.org/Tools/VEP
[2]: http://www.ncbi.nlm.nih.gov/clinvar/
Question: How do I read an mpileup file into Python?

I have tried several times, but it is just not working. I have a txt file instead on an mpileup file.

Answer: import pysam

# Read in a BAM file
bamfile = pysam.AlignmentFile("outputFile.bam","rb")
for line in bamfile:
print(line)

# Read in a SAM file
samfile = pysam.AlignmentFile("outputFile.sam","r")
for line in samfile:
print(line)

#Read in a CRAM file.
CRAMfile = pysam.AlignmentFile("outputfile.cram", "rc")

#Create the iterable and print the first line
it = iter(bamfile)
print(next(it))

#itt = iter(samfile)
print(next(itt))

# A different way to open the mpileup file.
f = open("/Users/m.o.l.s/outputFile.mpileup","rt")

line_1 = f.readline()
print(line_1)

# Split up the strings based on tabs
individuals = line_1.split("\t")
print(individuals)

# The value we are interested in is in the 3rd index position
The_Coverage = individuals[3]
print(The_Coverage)

# If this number is less than 10, the row needs to be displayed
if The_Coverage <= "10":
print(The_Coverage)
else:
print("The Coverage is not less than or equal to 10")
Question: Hello all. I have been having issues on running a python script. I have tried to resolve it by uninstalling and reinstalling pysam since the error is about pysam mostly. below is the command and error i get.

$ python2.7 overlapping_reads.py
Traceback (most recent call last):
File "overlapping_reads.py", line 4, in
import pysam
ImportError: No module named pysam

Can someone help?also, belo is portion of my script if that helps. Thanks

import argparse
from collections import defaultdict

import pysam

def Parser():
the_parser = argparse.ArgumentParser()
the_parser.add_argument(
'--input', action="store", type=str, help="bam alignment file")
the_parser.add_argument(
'--minquery', type=int,
help="Minimum readsize of query reads (nt) - must be an integer")
the_parser.add_argument(
'--maxquery', type=int,
help="Maximum readsize of query reads (nt) - must be an integer")
the_parser.add_argument(
'--mintarget', type=int,

Answer: > I have tried to resolve it by uninstalling and reinstalling pysam

Make sure you install pysam for the right python version. In your case, I suggest:

python2.7 -m pip install pysam
Question: I am wondering if the number of digits in NCBI Taxonomy ID have any meaningful concept? As I checked the 3 and 4 digits refer to group and species and more than 4 digit Ids refer to the strains? Is that true?

For example:

Tax ID: 536 [**Chromobacterium violaceum**][1]

Tax ID: 243365 Chromobacterium violaceum ATCC 12472

[1]: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=536&lvl=3&lin=f&keep=1&srchmode=1&unlock

Answer: I don't think so. Prokaryote taxa for which pure cultures have existed for a long time can be expected to have relatively few digits in their IDs..

Are taxids stable?

Taxids are assigned sequentially as taxa are added to the taxonomy database, and are never reused for a different taxon once they have been assigned. When two taxa are determined to be synonymous they are merged, and the taxid that disappears becomes a ‘secondary’ taxid for the one that remains. Taxa are occasionally removed from the taxonomy database (particularly internal nodes, during a taxonomic revision) – these taxids are deleted, and are not reused. Lists of merged and deleted taxids are included in the taxonomy dump files on the ftp site.

The public web pages at the NCBI only show taxa that are linked to public sequence entries. If a taxid is linked to a public sequence entry that is later removed, or assigned to a different taxon, the corresponding taxid may disappear from the NCBI web pages without appearing in either the merged or deleted taxid list.

http://www.ncbi.nlm.nih.gov/books/NBK54428/
Question: My multiple fasta looks something like that:

>chr2_1000-1020

NNNNNNNNNNNNNNNN

>chr2_1-20

NNNNNNNNNNNN

>chr2_20-40

...........

>chr2-200-220

And I want to sort all the sequences in numerical order. Any suggestions?

Answer: The following is not the best answer, but I cannot help mentioning my alpha-numeric sort, which I use daily. It is modified from unix sort and sorts "chr10,chr2" to "chr2,chr10" and also coordinates in the right way. You can use Pierre's command line to turn fasta into one line, I will just use bioawk for simplicity:

bioawk -c fastx '{print}' in.fa | sort -k1,1N | awk '{print ">"$1;print $2}'

where "sort" here is compiled from by modified sort. For me, this alpha-numeric sort has reduced a lot of tedious works related to coordinate sorting.

----------

Edit 2.3.21 by https://www.biostars.org/u/25721/

In the original answer it was `awk -c fastx` instead of `bioawk`, I made that change now.

One might also try `sort -k1,1V` to get a simple natural sort here in case that is of interest.
Question: Why do I get different number of lines when I use bash sort vs Bedtools sortBed

$ cut -f 1-4 BED.bed | wc -l

295658

-

$ cut -f 1-4 BED.bed | sort -k1,1 -k2,2n | uniq | wc -l

236207
_

$ cut -f 1-4 BED.bed | sort | uniq | wc -l

236207
_

$ cut -f 1-4 BED.bed | sortBed | uniq | wc -l

289619

-

$ cut -f 1-4 BED.bed | head -n 5

chr1	703951	761450	DUP
chr1	839849	839969	DEL
chr1	839858	839956	DEL
chr1	839859	839919	DEL
chr1	839872	839933	DEL

Answer: By default, bedtools sorts by chr and start, ignoring the end-coordinate. Therefore, it might happen that you get something like this:

chr1 4 15

chr1 4 17

chr1 4 15

Uniq iterates from the top to the bottom of this file. Therefore, it recognizes duplicated reads only, if they appear right one below the other, which is not the case in the above example.

I recommend you not to use bedtools for sorting, as indicated on the bedtools::sort manual page. Unix sort is faster and more memory efficient. If you want to combine sorting of your bed with deduplication, you can use the following command:

sort -k1,1 -k2,2n -k3,3n -k6,6 -u input.bed > output.bed

This command takes chr, start, end and strand into account, which are the essential information to describe a unique fragment. The -u then acts on all the columns, that were provided in the command, generating a deduplicated file.

If you are on a Mac, you may want to get the GNU core utilities e.g. via homebrew, as GNU sort provides a nice --parallel option for multi-threading.
Question: Hi Biostars,

I am looking for a tool for finding CpG islands in a draft genome-assembly. I see some tools which are web-based with input limit. I also see some tools that are specialised for mammalian genome promoters. The organism is non-mammalian. Some tools use epigenetic marks. I don't have DNA methylation data.

Is there a tool or a script that can use large `FASTA` files and find CpG islands?

Answer: Hi Satya,

Have you already tried [cpgplot][1], available [here][2]?

[1]: http://emboss.open-bio.org/rel/rel6/apps/cpgplot.html
[2]: ftp://emboss.open-bio.org/pub/EMBOSS/emboss-latest.tar.gz
Question: Hello all,

When working with publicly available data one often has to first find them in GEO, and then download the raw reads from the SRA archive database. I've been using the Entrez Direct tools from NCBI for this.

For example, using Entrez Direct tools from NCBI, I can put in the following query:

esearch -db sra -query GSM1467783 | efetch -format runinfo

And get the following output:

Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR1539207,2015-07-22,2015-12-01,37427735,1871386750,0,50,1424,,https://sra-download.ncbi.nlm.nih.gov/srapub/SRR1539207,SRX672143,,RNA-Seq,cDNA,TRANSCRIPTOMIC,SINGLE,0,0,ILLUMINA,Illumina HiSeq 2000,SRP045352,PRJNA257777,2,257777,SRS676636,SAMN02978909,simple,9606,Homo sapiens,GSM1467783,,,,,,,no,,,,,GEO,SRA178460,,public,B5573983CBB9C5E2046EB16D1C15A72F,25EC9BD6214D8EC0E7B270D41E6E861F

Now, I'm looking for a way to do it offline - to generate a big reference of GSM-to-.sra file correspondence, so by specifying the GSM ID or several, I can download all the relevant .sra file from the ftp.

I've found some R package wrappers around SQL databases; however, I'm puzzled that no table in SRAdb (I believe I've checked them all) includes GSM ids. I also haven't found the SRX/SRR ids in the GEOmetadb tables, although I might have missed something there.

At any rate, is there any database that establishes the correspondence between GSM and SRX/SRR identifiers?

If anybody can help me out with this, I'd be very grateful.

Answer: Have a look at the [SRA ftp server](ftp://ftp.ncbi.nlm.nih.gov/sra/reports/Metadata).
They provide a file called [SRA_Accessions.tab](ftp://ftp.ncbi.nlm.nih.gov/sra/reports/Metadata/SRA_Accessions.tab) which links various identifiers to each other.

With

```
grep ^SRR SRA_Accessions.tab | grep GSM
```

you can extract the SRR to GSM mapping.
Question: Hi friends,
I have some significantly differentially genes that are down regulated(negative fold change).
In the pathway analysis for these genes, is it correct interpretation to say that pathways related to these genes are downregulated in my samples as these genes are downregulated and in normal case if they were not down regulated they would be able to enrich the related significant pathways)

Answer: It only means that you might have something going on for that pathway in your data. What exactly is on you to figure out biologically. Say a pathway had 100 genes, and you find 10 of them being upregulated, but the molecular function of those 10 genes was to actually suppress (so negatively regulate) the pathway, then the pathway biologically would go down, even though these genes go up. Be careful with these enrichment tools, all they do is statistical enrichment analysis, the cell itself does not care about statistics. After DEG it requires biological knowledge of your system and downstream experiments to figure out the actual in vivo process. After all DEGs are just transcriptional, there is no guarantee this translates to actual protein activity changes in a pathway. Keep all that in mind.
Question: Hi everyone,

I'm trying to run Mutect2 for WES cancer data.

However, since their Resource bundle only supports h19 seems I cannot proceed (I want to compare it with Strelka2 results).

I've been looking for some hg38 interval_list file and I found:
`hg38_v0_HybSelOligos_whole_exome_illumina_coding_v1_whole_exome_illumina_coding_v1.Homo_sapiens_assembly38.targets.interval_list`

However, when I run the GenomicsDBImport I get the error (no matter if I use my own hg38 reference and .dict or the ones from GATK Resource Bundle):

A USER ERROR has occurred: Badly formed genome unclippedLoc: Contig chr1 given as location, but this contig isn't present in the Fasta sequence dictionary

So, my questions are:

1. Does anyone know if there is any release date for this hg38 based exome interval file?
2. Or the file I put is ok and the error is coming from somewhere else?

Answer: For the ones interested:

This error was a product of incorrect indexing of Fasta reference.

I used the hg38 reference available on the NCBI web.

However, all the chromosome regions are listed as `NC_10000.0` codes instead of 1, 2, 3 or chr1 chr2 and chr3.

I had to change it to another reference (I downloaded both .FA and .FAI from GATK Resource Bundle but it seems to work fine with hg38 from Ensembl too).
Question: Hi,

In the samtools manual, for the DP4 field it states that the:

"...sum can be smaller than DP because low-quality bases are not counted."

I wonder, what is the cut-off for this? How low does the quality of a base need to be and what determines whether a base is low-quality?

Thanks,

Jon

Answer: This is the mpileup -Q option, which by default is 13.