cannin/biostars_qa huggingface dataset:
Question: When I use snakemake, I try to calculate md5 value for several files. The files are as below:
```
$ ls -l|cut -d" " -f1,2,5-
-rw-r--r-- 1 0 Oct 21 19:46 A.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 B.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 C.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 D.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 E.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 F.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 G.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 H.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 I.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 J.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 K.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 L.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 M.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 N.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 O.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 P.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 Q.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 R.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 S.fastq.gz
-rw-r--r-- 1 131 Oct 21 19:46 test.py
-rw-r--r-- 1 0 Oct 21 19:46 T.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 U.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 V.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 W.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 X.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 Y.fastq.gz
-rw-r--r-- 1 0 Oct 21 19:46 Z.fastq.gz
```

So I created a Snakefile file called `test.py`:
```
rule md5:
input:
"{sample}.fastq.gz"
output:
"{sample}.md5"
shell:
"md5sum {input} > {output}"
```

and then I run the code and got the error as below:
```
$ snakemake -c -s test.py *gz
Building DAG of jobs...
MissingRuleException:
No rule to produce I.fastq.gz (if you use input functions make sure that they don't raise unexpected exceptions).
```

I don't know why this error occurs and how to solve this.

I would be appreciated if anyone could help me figure out it.

Thanks a lot.

Answer: Hello, you could pass the output instead of the input and let snakemake figure it out where the fastq file is located:

# -n is for a dryrun and -p is for snake to print the command line
snakemake -c -s test.py A.md5 -np

But this will work only for the A.fastq.gz, to run all you could do some shell string manipulation:

snakemake -c -s test.py $(ls *gz | sed -r 's:\..+::;s:$:.md5:')

This works but it is not very flexible, use glob_wildcards and a rule all.
Edit the test.py file:

from pathlib import Path

sample_dir = Path('.') # local dir

SAMPLES, = glob_wildcards(sample_dir / '{sample}.fastq.gz') # Extract the name of the samples to create wildcards.

rule all:
input:
md5s_files=expand("{sample}.md5", sample=SAMPLES) # force the creation of the files

rule md5:
input:
"{sample}.fastq.gz"
output:
"{sample}.md5"
shell:
"echo md5sum {input} > {output}"

then run:

snakemake -c -s test.py -np

You can further improve your code by passing the samples_dir variable as an argument inside a config.yaml file or through the command line.
Question: I have downloaded the NCBI nt database using the `blastdb_update.pl` perl script, but I want to blast some query file not on the whole nt database but on specific species. I know that when using blast locally it is possible to subset the nt/nr database using a list of GI identifiers, as explained [here][1].

However, [NCBI is phasing out GIs][2] and we should instead use accession.version identifiers. I have downloaded those for my species, below is part of the file `mygi.txt`.

When I run

`blastdb_aliastool -gilist mygi.txt -db nt -out sthg.out -title sometitle`

I obviously get

`BLAST Database error: Specified file is not a valid GI/TI list.` since I am not providing a GI list.

I cannot find any command-line option in the manual to specify that I want to filter the nt database by accession number; any idea of how I can achieve that? I bet this option will have to be added by the BLAST team at some point :)

-----------------------
`mygi.txt` below

AF324813.1
AF324814.1
AF324815.1
AF324816.1
AF324817.1
AF324818.1
AF324819.1
AF324820.1
AF324821.1
AF324822.1
AF324823.1
AF324824.1
AF370451.1
AY198341.1
AY198342.1

[1]: https://www.biostars.org/p/6528/#6532%29
[2]: http://www.ncbi.nlm.nih.gov/news/03-02-2016-phase-out-of-GI-numbers/

Answer: This solution adds a step but until NCBI updates the `blastdb_aliastool` to accept accession numbers this may the only way.

You can use `blastdbcmd` from blast+ package to retrieve sequences from `nt` db as fasta file followed by `makeblastdb` to make the blast indexes for the subset of sequences. my_acc.txt file is the file with accession numbers (one per line).

blastdbcmd -db /path_to/nt -entry_batch my_acc.txt -out my_seq.fa
Question: I have a bunch of bedgraph files with very little annotation information, and I have no idea whether the data was generated from whole genome or targeted bisulfite sequencing. What is the best way to determine this?

Answer: Open the files in IGV and check if the reads align to the whole genome or to only a subset of the genome (e.g. the exome)?
Question: I am trying to gain insight into the potential function of a list of candidate effector protein sequences I have identified. I have searched using NCBI BLASTP against the nr database, performed GO analysis, searched for known domains using InterProScan and blasted them against the PHI-Base database, but still, the majority are annotated as “hypothetical proteins” and have no known protein domains.

I am struggling to decide where to go next to try to gain insight into their potential function. I have considered trying to predict the 3d structure, but I would have to do this for 72 sequences and then what? I can’t find a database to search protein structures against.

Does anyone have any suggestions where to go next? Are there other papers, tools, or approaches you could suggest to gain more insight into potential function?

Thanks,
Jamie

Answer: Though unlikely, it is possible that all your proteins are uncharacterized and have no known domains. It is possible that some of them are not proteins at all. You may want to try the links below if you haven't already done so.

- https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi
- https://pfam.xfam.org/
- http://smart.embl-heidelberg.de/
- https://toolkit.tuebingen.mpg.de/tools/hhpred
Question: So I have managed to successfully query 1000 Genome VCF files before by following tutorials, I thought I could apply this to the COSMIC database. I can't.

I'm trying to make use of some WGS data from COSMIC. As I don't know the file structure I ran the following in an attempt to view the headers:

tabix -H ftp://ngs.sanger.ac.uk/production/cosmic/wgs/CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz
[get_local_version] downloading the index file...
[kftp_connect_file] 550 No such file.
[download_from_remote] fail to open remote file.
[tabix] failed to load the index file.

...so then I downloaded the file:

tabix -H CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz
[tabix] the index file either does not exist or is older than the vcf file. Please reindex.

How do I reindex a vcf file? I tried this (I don't know what I'm doing at all, sorry):

tabix -f CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz
[tabix] was bgzip used to compress this file? CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz

....and so in a last ditch attempt I tried just running vcftools (I thought this would filter out any mutations not on X):

vcftools --gzvcf CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz --chr X --out test

VCF index is older than VCF file. Will regenerate.
Building new index file.

Reading Index file.
File contains 641910 entries and 0 individuals.
Filtering by chromosome.
(list of chromosomes)
Skipping Remainder.
Keeping 28481 entries on specified chromosomes.
Applying Required Filters.
After filtering, kept 0 out of 0 Individuals
After filtering, kept 28481 out of a possible 28481 Sites

So didn't work as expected either. I then tried running tabix again but it still gives me the error that an index file doesn't exist.

Are there any read me files or guides that I've just not been able to find? I literally have no idea what I'm doing.

Answer: Roughly, you'll need to do these steps (untested):

gunzip CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz
bgzip CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf
tabix -p vcf CosmicCodingMuts_v64_26032013_noLimit_wgs.vcf.gz
Question: HI,

I got a question about how to assign color value in image() function. I thought the col option works as a vector with the same length of Z, with color of each cell in the matrix Z corresponding to the values in. However, it confused me with different result.

To make my question simple, let's say I want to make a color image for a vector `x=1:4` with corresponding color for each cell.

```r
x=1:4
image(1,1:length(x), matrix(x, nrow=1, ncol=length(x)), col=c("blue","red",'green','yellow'))
```

The output image looks like this:

![output image 1][1]

However, if I use a different x, say

```r
x=c(3,1,2,1)
image(1,1:length(x), matrix(x, nrow=1, ncol=length(x)), col=c("blue","red",'green','yellow'))
```

I expected the color (from bottom to up) order should be green-->blue-->red-->blue, however, I got a very different order:

![output image 2][2]

I guess I misunderstand how to use the color option in the function. Can someone please help to advise? Thanks

-Xianjun

[1]: https://www.dropbox.com/s/vzzt4t8u0nayy85/image1.png?raw=1
[2]: https://www.dropbox.com/s/g0traftglvtknrg/image2.png?raw=1

Answer: OK I figured it out myself:

Unless the breaks is set, the range of Z (or x in above example) is evenly cut into N intervals (where N = the length of color) and values in Z are proportionally assigned to the nearest color.

For example, when `x=c(3,1,2,1)` and `col=c("blue","red",'green','yellow')`, the minimal of x is assigned as the first color, and max to the last color. Any value between is calculated proportionally to a color. In this case, 2 is the the middle one, according to the principal that intervals are closed on the right and open on the left, it's assigned to "red". So, that's why we see the colors are yellow-->blue-->red-->blue.
Question: Hello,

I want to compile a BED file with all genomic regions where I find a given sequence (say, I have a 10-nucleotide sequence, and 10^9 nucleotides in the genome). I have previously tried solving this task with a weight matrix based approach, but it takes too long. So I decided to simplify the task: just a fixed sequence, no mismatches. I am thinking of something like Bowtie, but with all non-unique hits reported. Looks like a pretty simple task. Did someone tried solving this?

Thanks!

Answer: It'd probably be faster to just use Biopython. Iterate over the entries in a fasta file, and use find iteratively to produce all of the found regions.
Question: Hi all,
I'm trying to use Diffbind to identify differentially open regions (identified using ATAC-seq) between three different cell types. I have three biological replicates for each cell type.
I'm building now the sample sheet to import the data into Diffbind, but I don't know what I should indicate as a control since ATAC do not have controls. Is there a way to avoid this problem?

Answer: Control reads, such as those from an Input sample, are optional in DiffBind, so you can just leave out the control column in the sample sheet.
Question: I have an alignment of an assembled contig with the reference genome. The contigs are very long (PacBio + HiC assembly), and one "query" will generate multiple alignments with the "reference", as a consequence of structural variation, which breaks up the alignments.

I want to figure out if those alignments are in *the right order* relative to the contig they originate from. Next step would be to check if any of these alignments have deletions - leading to gaps, or change direction (due to inversions). I can't find any pysam attributes which seem do to the trick, right now, and I'm not sure if it's encoded in some obscure tag of the bam?

Suggestions?

Answer: Hi,

We (in our lab) had a similar problem, were we have reads with many supplementary alignments. Apparently, more people has the same problem as [you and me][1]

In order to get the supplementary tag as a function of the query (read), what I did was the following:

1. Align the reads with minimap2 (I have nanopore data), and output in paf format
2. Sort the data by column 1 and then by column 2 (sort -k 1,1 -k 2,2). In this case, you are sorting the alignments first by query name (you get all the alignments for a read one after each other) and then when you sort by column two you get the supplementary alignments one after each other **as a function of the query**

Let me know if this helps,

[1]: https://github.com/lh3/minimap2/issues/386
Question: Hi,

I'm looking for official algorithm documentation on Picard MarkDuplicates and SamTools rmdup, but I can't find it. I have found numerous posts in "google land" where people state why one is better, but I want to know exactly how they both work (preferably without going through the code). For example, I have "heard" that MarkDuplicates is more "intelligent" because it allegedly considers variants within a read rather than just looking at where reads begin.

Can anyone point me towards a paper or documentation that discusses the true differences between the algorithms?

Thanks for your help!

Mark

Answer: warning this answer is old. check the new version of samtools.

SamTools rmdup 'only' compares two reads on chrom and pos (which could be wrong if two reads come from two different libraries) and **removes** reads from the BAM: information is lost.

picard sets the sam flag 1024 but does not delete the reads. two pairs of reads are compared , as far as I know, using the chrom, the pos, the group-id (sample...) + (flowcell , lane, X,Y for optical dups) (,and the cigar string ?).
Question: I have a huge dataset of SNPs that I am trying to get hg19 locations for based only on the rsID. Right now I am doing that by just downloading the latest version of dbSNP, turning it into an sqlite database, and doing a huge long running query. This works, although it takes forever, but it has the problem that it fails for all SNPs that have been merged, e.g. rs111199278.

I need some way in batch to be able to get the position of these old rsIDs, possibly by just converting them to the latest rsID in dbSNP and then looking them up with my slow database query.

Is there a good way to do this efficiently? I think I will end up needing to lookup about 100K-500K rsIDs (I don't know yet because my big location lookup hasn't finished).

Answer: download and insert the content of ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/database/organism_data/RsMergeArch.bcp.gz into a sqlite database ?

$ curl -s "ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/database/organism_data/RsMergeArch.bcp.gz" | gunzip -c | grep -m1 -w 111199278
111199278	3095314	144	1	2015-04-01 13:45:00.837	2015-07-14 23:14:01.543	3095314	1	rsm

PS: I don't think sqlite3 is a database of choice for such large dataset.
Question: I have a very beta but working version of my software that might be a year away from full completion and two years away from completing analyses with the software that solves a novel biological problem (hence the paper is ~2+ years away; I do not just want an application note).

If I want to create a Python package and release the code on my github, what are the risks and benefits? I can see several benefits: users, feedback, pull-requests and bugfixes. But what are the downsides? Can it hurt me in any way?

For maximum awkward, the software is just a reimplementation of an old, but widely used algorithm with a software package that is unlicensed (as far as I can see). See here for more info: https://www.biostars.org/p/136001/

**Edit**: is a license enough to protect me (force others to give me rightful credit)? From a protection viewpoint, is there no benefit in writing and publishing a preprint? I'd rather avoid the latter for now; a README.md is good enough for me.

Answer: There are more pros than cons. Sharing your beta code and allowing other users to play with it give you feedback, free beta testers and, luckily, new ideas and suggestions (an extra point of view is always fresh and enriching). Additionally, sharing your code could bring you some publication references talking about (free advertising!) and increasing the number of people knowing it.

About the downsides, GitHub allows you to set licenses for your code (MIT License for example), with this you can avoid the only drawback I have found: somebody using your software without referencing it. The license includes the date so your software will always be "protected" against posterior inappropriate use.

About being a reimplementation, there will be no problem while you will always reference the original author.
Question: Hello all,

Is it fine if I stitch all the 167 contigs from my genome into 1 single fasts file?

Are there any tools which can be employed for this purpose.

I have a strain of bacteria with 167 contigs and I want to stitch them up together under 1 single FASTA header.

Thanks a lot

Cheers
Optimist

Answer: It's not illegal, but that was possibly not what you meant ;) For most other use-cases it is definitely **not** fine. I'll just give you some hints:

- You cannot tell which sequence is belonging to which contig anymore
- You don't know the natural order of contigs
- You are creating artificial chimeric sequences at the transitions
- Sequences can come from different replicons
- ....

- You can possibly try to scaffold the contigs if there are additional long-range sequencing libraries (Fosmid, pacbio, ...) (SSPACE is the program)
- You can keep the contigs as they are in a single multi-fasta file
- If you still think you need to do this, you need to re-think the reason. There is probably a problem with your approach.
Question: I am trying to identify runs of homopolymers in sequences of a bed and have a perl script that works for the desired output below. There are ~55,000 lines in the custom.bed, is there a way to download all the sequences in the bed in the desired output format? UCSC table browser will give the output in the correct format however there is a limit of 1,000 entries per query. The below seems close but the output is formatted incorrect. Thank you :).

custom.bed
chr1	948953	948956	chr1:948953-948956	.	ISG15

BED=custom.bed
for chr in `seq 1 22` X Y
do
wget -O - -q http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr1.fa.gz | gunzip -c >> hg19.fa
done

fastaFromBed -fi hg19.fa -bed $BED -fo $BED.fasta

current output:
>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

desired output of each line for fasta sequence:

>hg19_refGene_NM_002335_0 range=chr11:68080173-68080283 5'pad=10 3'pad=10 strand=+ repeatMasking=none
gccggacaacATGGAGGCAGCGCCGCCCGGGCCGCCGTGGCCGCTGCTGC
TGCTGCTGCTGCTGCTGCTGGCGCTGTGCGGCTGCCCGGCCCCCGCCGCG
Ggtaggtgggc

Answer: Hello, cmccabe.

Thank you for your question about getting FASTA sequences from the UCSC Genome Browser.

Are you uploading your "custom.bed" as a custom track and then attempting to get the FASTA sequence from the Table Browser? You are correct that there is a limit of 1000 regions if you are using the Table Browser's "define regions" option, but for custom tracks, there is no such limit.

Here are some basic steps that you could use to get the sequence you're interested in:

1. Upload your "custom.bed" as a custom track: https://genome.ucsc.edu/cgi-bin/hgCustom?db=hg19
2. After uploading, navigate to the Table Browser:https://genome.ucsc.edu/cgi-bin/hgTables.
3. Make the following selections:
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Custom Tracks
track: My Custom Track
table: myCustomTrack
output: sequence
output file: enter a file name to save your results to a file, or leave blank to display results in your browser

4. Click "get output".
5. Select your sequence output options.
6. Click "get sequence".

You should end up with a FASTA with names similar to what you are looking for in your "desired output".

If you have any follow up questions, it would be helpful if you could post them to our Google Groups forum: https://groups.google.com/a/soe.ucsc.edu/forum/#!forum/genome, that way our whole team can see the question and help with an answer.

Thanks,

Matthew from the UCSC Genome Browser
Question: Hi.

I have a question from using pipeline that find the structural rearrangements like SVDetect or BreakDancer.

I ran above software and now I am having a results from SVDetect and BreakDancer.

(Input files is a tumor bam file based on paired end and output files is just a txt file.)

Results have many columns that describe the sequence information like start position, end position, class of

type(indel, translocation,duplication and etc...).

I almost understand why these types are classified by software.

However, I want to see manually that really my tumor bam file has the structural variation which is previously

classified by software.

For example, I used SVDetect which find structural variants and got a output result and its content is below.

chr10 255485 255896 chr10 255738 256047 (F,F,F,F,F,F) (R,R,R,R,R,R) TRANSLOC

As I mentioned, I already understood the information above conceptually. but my ultimate goal is just want to identify whether TRANSLOC really happened at above position in tumor bam files or not.

Because many software has mis classification. Therefore, It is required for me to see variants manually not only seeing the output file.

Anyway my question is how you validate the output files ?

If you have your own validate process.. could you explain me in detail?

I heard the IGV program and I already used this tool but I can't find anything information that i want.

I am looking forward to you reply.

Thank you

Answer: IGV is definitely a good tool for this task. There are some details that might be applicable [here][1]

[1]: http://www.broadinstitute.org/igv/interpreting_pair_orientations

Just FYI, running the structural variant detection on the "normal" sample from the same patient or a pool of "normals" can very greatly reduce your false positive results if you exclude any SV in the tumor that is overlapping one in the normal (or normals).
Question: Hello,

I'd like to remove a tag (not a read) from a BAM file. I would like to remove the OQ tag that keeps original qualities after running BQSR. Any suggestions on how I might do this?

Thanks,
Mark

Answer: sed 's/\tOQ\:Z\:[^\t]*//' in.sam > out.sam