cannin/biostars_qa huggingface dataset:
Question: I am using Tophat2 (v2.1.1) to find fusions on RNA-seq files.

To test the --fusion-search function, I launch the example in the manual :

tophat -o tophat_outpu2 -p 8 --fusion-search --keep-fasta-order --bowtie1 --no-coverage-search -r 0 --mate-std-dev 80 --max-intron-length 100000 --fusion-min-dist 100000 --fusion-anchor-length 13 --fusion-ignore-chromosomes chrM /GenomeRef/hg19 SRR064286_1.fastq SRR064286_2.fastq

And everytime I got the same error :

Error running /programs/tophat/tophat-2.1.1/install/bin/tophat_reports --min-anchor 8 --splice-mismatches 0 --min-report-intron 50 --max-report-intron 100000 --min-isoform-fraction 0.15 --output-dir tophat_output2/ --max-multihits 20 --max-seg-multihits 40 --segment-length 25 --segment-mismatches 2 --min-closure-exon 100 --min-closure-intron 50 --max-closure-intron 5000 --min-coverage-intron 50 --max-coverage-intron 20000 --min-segment-intron 50 --max-segment-intron 100000 --read-mismatches 2 --read-gap-length 2 --read-edit-dist 2 --read-realign-edit-dist 3
...
tophat_output2/insertions.bed tophat_output2/deletions.bed tophat_output2/fusions.out tophat_output2/tmp/accepted_hits tophat_output2/tmp/left_kept_reads.mapped.bam,tophat_output2/tmp/left_kept_reads.candidates tophat_output2/tmp/left_kept_reads.bam tophat_output2/tmp/right_kept_reads.mapped.bam,tophat_output2/tmp/right_kept_reads.candidates tophat_output2/tmp/right_kept_reads.bam

./SeqAn-1.4.2/seqan/basic/basic_exception.h:236 FAILED! (Uncaught exception of type St12out_of_range: basic_string::substr)

Dis someone can help me with this error ?

Thanks.

Answer: The right way to do is using Tophat2 v2.1.0...........
so simple but useful
Question: I've tried to search a bit first, but have been unable to find a quick answer. In the documentation for biopython I come across this:

#Distance is calculated based on the number of shared 6mers.

This is the default alignment algorithm for MAFFT. What is a 6mer?

Also when should I choose the alternatives over the 6mer pairwise alignment? Those alternatives are:

Needleman-Wunsch (global pairwise)

Smith-Waterman (local pairwise)

Local pairwise with generalized affine gap cost (Altschul 1998)

Answer: A 6-mer is a sub-sequence of 6 letters (nucleotides or amino-acids).

For example, a length 10 DNA sequence has 5 possible 6-mers: the sequences using bases 1..6, 2..7, 3..8, 4..9, and 5..10. Think of a sliding window of length 6 moving across the sequence one letter at a time.

Aligning two sequences with the alternatives you list is the "optimal" way but it takes time proportional to the product (*) of the lengths of the two sequences. The distance can be measured from the alignment score.

To speed this up for the case where the sequences are reasonably similar, another way is to *count the number of shared k-mers*. MAFFT must use k=6. This can be done more quickly than full alignment.

The decision is a trade-off between computation and sensitivity. k=6 might be ok for protein but may fail for DNA if the sequences diverge to any great extent.
Question: I'm using geneBody_coverage2.py in RseQC package to calculate Read coverage over gene body, which need a bed12 file. I did alignment based on Homo_sapiens.GRCh38.84.gtf , so i need a Homo_sapiens.GRCh38.84.bed. But gtfToBed can't produce a BED12 file, how can I get it? anyone can help me?

Answer: `awk '{if($3 != "gene") print $0}' Homo_sapiens.GRCh38.84.gtf | grep -v "^#" | gtfToGenePred /dev/stdin /dev/stdout | genePredToBed stdin Homo_sapiens.GRCh38.84.bed`

You can download the two tools from UCSC.

If you're curious, this is from [here](https://github.com/maxplanck-ie/data_repository/blob/master/organisms/Makefile.inc#L56-L59) and is essentially the code I use to generate the reference data for our institute.
Question: How can I set default values for cwl files?

Here is my `test.cwl` file

cwlVersion: v1.0
class: CommandLineTool
baseCommand: echo
inputs:
example_flag:
type: boolean
default: true
inputBinding:
position: 1
prefix: -f
example_string:
type: string
default: "abc"
inputBinding:
position: 3
prefix: --example-string
example_int:
type: int
default: 5
inputBinding:
position: 2
prefix: -i
separate: false
example_file:
type: File?
inputBinding:
prefix: --file=
separate: false
position: 4

outputs: []

With

$ dockstore tool convert cwl2json --cwl test.cwl

I get the following output

{
"example_string": "fill me in",
"example_file": {
"path": "fill me in",
"class": "File"
},
"example_flag": false,
"example_int": 0
}

Where are the default values I set in the `test.cwl` file?

Answer: As a heads-up, this is now a feature for Dockstore 1.1. You should be able to do the following

$ dockstore tool convert cwl2json --cwl test.cwl
{
"example_string": "abc",
"example_file": {
"path": "fill me in",
"class": "File"
},
"example_flag": true,
"example_int": 5
}
$ dockstore tool convert cwl2yaml --cwl test.cwl
example_string: abc
example_file:
path: fill me in
class: File
example_flag: true
example_int: 5

$ dockstore tool convert cwl2yaml --cwl test.cwl > test.yaml
$ vim test.yaml
$ dockstore tool launch --entry test.cwl --local-entry --yaml test.yaml
Question: How I can replace number in R 1 to11, 2 to 22 and 0 to 10, also -9 as NA?

id rs143 rs148 rs149 rs1490
1 02003s -9 1 2 1
2 02003s -9 0 1 2
3 02003s -9 1 1 2
4 02003s -9 0 1 1
5 02003s -9 0 1 1

I want to replace numbers as this format:

id rs143 rs148 rs149 rs1490
1 02003s NA 11 22 11
2 02003s NA 10 11 22
3 02003s NA 11 11 22
4 02003s NA 10 11 11
5 02003s NA 10 11 11

I tried this

df[df==1]<-11
df[df==2]<-22
df[df==0]<-10

But this error occur:

In [<-.factor`(`*tmp*`, thisvar, value = 22) :
invalid factor level, NA generated

Answer: As mentioned in the comments some of the columns are class of *factor*, not *numeric*. If all SNP columns are numeric, then your existing code should work.

Try below code to convert, factor **and** numeric columns into numeric:

# lookup for allele recoding
lookup <- setNames(c(NA, 10, 11, 22), c(-9, 0, 1, 2))

res <- cbind(df[, 1, drop = FALSE],
data.frame(lapply(df[, -1], function(i) lookup[ as.character(i) ])))

res
# id rs143 rs148 rs149 rs1490
# 1 02003s NA 11 22 11
# 2 02003s NA 10 11 22
# 3 02003s NA 11 11 22
# 4 02003s NA 10 11 11
# 5 02003s NA 10 11 11
Question: Hi Biostars community,

I have bam files with RNA-seq results from ~1000 of single cell samples. For every transcript in each sample, I would like to calculate what percentage of spliced transcript's length is covered by aligned reads. Besides, It would be great to have this value for coding part of transcripts only, or, ideally, for every exon.

What would be the most straightforward way to do this?

Thanks a lot!

Answer: Thanks to WouterDeCoster and Bastien Herve for their answers!

I have found out that bedtools does exactly what I need:

bedtools coverage -hist -abam "your bam file" -b "your bed file" > "file with results"
Question: Hi everyone,

I have the fastq files for some miRNA libraries prepared with the QIAseq miRNA Library Kit. I have to do the UMI extraction, but the problem is that the UMI is after a common sequence for all the reads, such as this:

> NNNNNNNNNNNNNNNNNNN**AACTGTAGGCACCATCAAT***XXXXXXXXXXXX*NNNNNNNNN

Where the N are the miRNA sequences, the **bold** part is the common sequence for all the reads and the part with all the X is the part with the UMI sequence.

How could I remove the bold part and append the UMI to the header of the fastq file?
The problem is that I have seen that around 3-5% of the reads don't have the common sequence, I suppose that there are sequencing errors and some part of this sequence is changed in some reads, but I don't know how to accept one letter change in the common part.

Thank you very much!

Answer: You should be able to do this with UMI tools using the regex UMI extraction mode.

you can do something like:

```
umi_tools extract --extract-method=regex \
--bc-pattern=".+(?PAACTGTAGGCACCATCAAT){s<=2}(?P.{12}).+" \
-I input.fastq.gz \
-S processed.fastq.gz
```
The `{s<=2}` means "allow two mismatches in the common sequence". Note that this will leave both the Ns at the start of the sequence and the Ns at the end of the sequence intact. If you wish to remove the Ns at the end then you the regex:
`.+(?PAACTGTAGGCACCATCAAT){s<=2}(?P.{12})(?P.{9})`

See more details here:
https://umi-tools.readthedocs.io/en/latest/regex.html#regex-regular-expression-mode
Question: I have a set of 500 100-base sequences. I want to make new sequences by randomizing this set into 500 new sequences such that the new set has the same dinucleotide distribution. So, for each sequence in my original set, I'm cutting the sequence up into mononucleotides, putting them in a hat, and drawing them out one by one, but the randomized result should have the same dinucleotide distribution. How do I do this?

The best idea I can come up with at present is to simply repeat a purely random shuffling until I find a result that is sufficiently close to the original. (closeness may be assessed with a paired t-test or, more simply, by taking the difference between the original distribution and the random one.) Is there a better way? I could slice the original sequences by twos to jump-start the process?

Background: I'm trying to reproduce a procedure described in in Weirauch et al. Nature 2012, (http://www.nature.com/nbt/journal/v31/n2/full/nbt.2486.html?WT.ec_id=NBT-201302). In trying to evaluate a TF binding prediction algorithm's success at predicting in vivo binding sites, they start with some ChIP-Seq or ChIP-exo data. From these measurements, they pick 500 loci where the TF is credibly believed to bind, and then test the algorithm's ability to distinguish the true sites from 500 control sites. This is one of the methods they use to create control sites, and they describe it thusly: "500 randomly shuffled positive sequences, where dinucleotide frequencies were maintained." This is the whole description in the methods section, leaving a few things unclear.

1. Do I shuffle the sequences individually or as a set? I mean, do I (a) slice all the sequences in the set and put them in a hat, and then test the dinuc distribution of the new set or (b) shuffle each sequence individually, as described above?
2. How precisely is the dinucleotide distribution maintained?

Answer: Most people use this script ([**altschulEriksonDinuclShuffle.py**][1]) to shuffle sequences while maintaining dinucleotide frequency (It's even [code][2] that's present in the MEME suite). written originally by P. Clote, or a derivative of it. I found this copy online in the Bias-Away package by the Wasserman lab. It works by shuffling each sequence to make a permuted sequence with the original dinucleotide content. So the entire procedure is to shuffle each peak sequence and compare TF binding scores on the original sequences versus all scrambled sequences.

[1]: https://github.com/wassermanlab/BiasAway/blob/master/altschulEriksonDinuclShuffle.py
[2]: http://web.mit.edu/meme_v4.9.0/bin/fasta-dinucleotide-shuffle
Question: I am trying to subset the samples according to condition type from a particular GSE from Gene Expression Omnibus (GEO). This particular GSE has data from two platforms:

- Affymetrix Human Genome U133A Array
- Affymetrix Human Genome U133B Array

I wanted to take samples from both these platforms and proceed with normalisation and deducing expression values for further WGCNA analyses.

My question is can we do the normalisation of samples from both these samples together or should I do it separately and later merge them (As the probe IDs of the two platforms will be different) ?

Thanking you in advance!

Answer: This has been done previously in 2 studies (possibly many more):

- Secreted
frizzled-related protein 1 loss contributes to tumor phenotype of
clear cell renal cell carcinoma
- Pathway
signature and cellular differentiation in clear cell renal cell
carcinoma

Both of these utilised this GEO dataset: GSE6344

The methods that they used for merging, which involves the determination of a scaling factor on raw signal data, are here: ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE6nnn/GSE6344/suppl/GSE6344_MethodsDetails.PDF

This topic has also been discussed on both Biostars ( https://www.biostars.org/p/56657/ ) and Bioconductor a few times.

*[just briefing everyone]*

--------------------------------
--------------------------

**The key to your particular question is that you just want to run WGCNA on the data.** As WGCNA is fundamentally based on correlation, in this case, you just need to process each array experiment independently and then combine the final normalised datasets together. The merged dataset then becomes the input to WGCNA. Obviously only combine the datasets where there are common genes.

An optional transformation for you is to normalise the arrays independently, covert the normalised expression values to the Z-scale, and then merge everything together.

If you want to start conducting statistical analyses other than WGCNA, then it's a different story and that particular topic is hot and provides for debate across the forums.

Kevin
Question: Hi all,

I want to count the number of SNP for each chromosome in the raw VCF file. What is the best idea?

Best Regard

Mostafa

Answer: UPDATE 2021: if your VCF is indexed: `bcftools index -s indexed.vcf.gz`

grep -v "^#" in.vcf | cut -f 1 | sort | uniq -c
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: Very basic, but it helped me get started, although I was mainly just interested in importing tabular data and plotting:

Zuur, Alain et al: A Beginner's Guide to R
Question: Hi!

I have a pile of unindexed reads, including some PhiX, from a poor run on a NexSeq. Can anyone suggest some tools or techniques I can apply to this data just to try get a better idea of why the run went poorly?

Much appreciated.

Answer: You could try blasting a few reads to get an idea whether it's poor quality, PhiX,...
There will be more sophisticated methods for sure but it's a start :)
Question: I have copies of the same `fastq file` (let's call it `samp1_R1.fastq.gz`) in two different locations (e.g., `/foo/samp1_R1.fastq.gz` and `/bar/samp1_R1.fastq.gz`). These files are supposed to be identical (the one at `/foo` is a "working copy" while the one at `/bar` is a backup).

I happened to generate MD5 checksums from both these files today, and to my surprise, the hashes are different:

$ md5sum /foo/samp1_R1.fastq.gz
13a75f5e319fa772faa85beb04317718 /foo/samp1_R1.fastq.gz

$ md5sum /bar/samp1_R1.fastq.gz
e7a4f3c293361fa569a192f8a26a141d /bar/samp1_R1.fastq.gz

Does this truly mean that the files are now different? I cannot think of anything that could have modified one of these files (the "working copy" was only ever used as an input to standard RNA-seq tools). Could the mere act of uncompressing (and recompressing) the file update the hash like that? (That's the only thing I can recall doing to the file itself.)

Now I'm also wondering whether I should rerun analyses that were based off of this file.

Your inputs would be appreciated.

Answer: different compression level ? names changed ?

check with:

gunzip -c /foo/samp1_R1.fastq.gz | md5sum

gunzip -c /bar/samp1_R1.fastq.gz | md5sum

or/and

paste <( gunzip -c /foo/samp1_R1.fastq.gz) <( gunzip -c /bar/samp1_R1.fastq.gz) | awk -F '\t' '($1!=$2)'
Question: Hello

I'm very confused, can any one brief me what is left align and parsimonious mean with simple example ?

Thanks!

Answer: This paper [Unified representation of genetic variants][1] or this [wiki page][2] from the same authors explain normalization very nicely.

[1]: https://academic.oup.com/bioinformatics/article/31/13/2202/196142
[2]: https://genome.sph.umich.edu/wiki/Variant_Normalization
Question: I am trying to run qualimap on multiple bam file via command.

qualimap multi-bamqc -d bam_file.txt -gff target.interval_list

But it gives an error as follows:

qualimap multi-bamqc -d bam_file.txt -gff target.interval_list
Java memory size is set to 120000M
Launching application...

OpenJDK 64-Bit Server VM warning: ignoring option MaxPermSize=1024m; support was removed in 8.0
QualiMap v.2.2.1
Built on 2016-10-03 18:14

Selected tool: multi-bamqc
Running multi-sample BAM QC

Checking input paths
Failed to run multi-bamqc
java.lang.RuntimeException: The raw data doesn't exist for sample: 001X1.sorted.deduped.bqsr.bam
Folder path:/media/molmed/Data/PJRCT1/bam
Please check raw data directory is present.

bam_file.txt conatins

`head bam_file.txt`

001X1.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
001X2.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
001M1.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
001M2.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
003X3.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam
003M1.sorted.deduped.bqsr.bam	/media/molmed/Data/PJRCT1/bam

I also tried by giving the complete path to the file as `/media/molmed/Data/PJRCT1/bam/001X1.sorted.deduped.bqsr.bam`

and also by adding `/` to the folder ends as `/media/molmed/Data/PJRCT1/bam/`. But all gives the same error

Answer: Summary:

To run `multi-bamqc` in `qualimap`, you should have a `bed` file containg six columns: You can see the comments [here][1]

cat bed3.bed | perl -lane 'print "$F[0]\t$F[1]\t$F[2]\t.\t0\t."' > target.bed #this will contain 6 columns

The file containing name of the sample and paths should be as:

sample_1 /media/User/Data//bam/sample_1.bam #file name with extension `.bam`
and run the following command

qualimap multi-bamqc -r -d bam_file.txt -gff target.bed

[1]: https://www.biostars.org/p/14889/
Question: I called structure variant using Sniffles for two samples (using PacBio reads), I would like to compare the resulted two files to know what is common and different in both files.

what first come to my mind is to use vcf-isec, but I do not know how it will deal with the deletion size if they are exist in the two files but different in length for example, translocation also if its origin is the same but the location transferred to is different! etc ..

vcf-compare could be used but I do not know if it is the right way!

is there is a tool that could do this in an efficient way? anyone have experience with that?

Thanks,

Answer: One simple way to do this is to merge the SV calls and look at which calls collapse. I've writing some code to do this in VCF format.

https://github.com/zeeev/mergeSVcallers
Question: Hello all,

I have a question about how to pass cnvkit's output to pureCN to account for the tumor cellularity,

According to https://bioconductor.org/packages/release/bioc/vignettes/PureCN/inst/doc/Quick.pdf, Page 6

I first do:

Rscript $PURECN/NormalDB.R --outdir $OUT_REF --normal_panel $NORMAL_PANEL \
--assay agilent_v6 --genome hg19 --force

When I run this, I am asked for --coveragefiles. Can I provide a file from cnvkit for this?

Secondly:

> cnvkit.py export seg $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.cns
> --enumerate-chroms \
> -o $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.seg

And finally:

> Rscript $PURECN/PureCN.R --out $OUT/$SAMPLEID \
> --sampleid $SAMPLEID \
> --tumor $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.cnr \
> --segfile $OUT/$SAMPLEID/${SAMPLEID}_cnvkit.seg \
> --normal_panel $OUT_REF/mapping_bias_agilent_v6_hg19.rds \
> --vcf ${SAMPLEID}_mutect.vcf \
> --statsfile ${SAMPLEID}_mutect_stats.txt \
> --snpblacklist hg19_simpleRepeats.bed \
> --genome hg19 \
> --funsegmentation none \
> --force --postoptimize --seed 123

Here, could someone please advise me on what snpblacklist is and what is being used as --normal_panel?

I tried reading the manuals but I am still confused.

Any suggestions would be appreciated, thanks!

Answer: Hi,

if your VCF file contains a matched normal, then --normal_panel is not crucial. Essentially, PureCN assumes that heterozygous SNPs have an expected allelic fraction of 0.5 without copy number events. By providing a pool of normals, PureCN can check all common SNPs for deviation from 0.5 and can adjust or discard poor quality SNPs. When a matched normal is available, PureCN has at least access to a single normal.

Have a look at the FAQ section how to create a panel of normal VCF. Then there is a NormalDB.R script that can generate the mapping_bias.rds file for you. But again, for test runs, you won't need this.

--snpblacklist can be any BED file. It will simply ignore every variant in the VCF that falls into those regions. Have a look at the main vignette where it downloads the simple_repeats track from UCSC.

And yes, PureCN accepts CNVkit coverage files.

If you use most recent GitHub (which will become new stable next week), you can replace --funsegmentation none with --funsegmentation Hclust. This might undo some over-segmentation when present.

Best,
Markus
Question: Hi everybody,
I have the counts (obtained by HTSeq) for a lot of genes(~58,000) at different time points (replicates).

gene t1_S1 t1_S2
ENSG00000000003.14 0 0
ENSG00000000005.5 0 0
ENSG00000000419.12 1 3
[...]

I woul like to calculate the correlation between the counts over each gene at the same timepoint to understand how reproducible the replication timing and progression is for each repeat.
Any suggestions?

Answer: Do you want to test the correlation between the different timepoints or between the different genes.

Let say you have 10 timepoints and 58000 genes

To test the different timepoints :

cor(xx, method="pearson")
will give you a 10x10 matrix , so 100 correlations calculation (even though I guess the `cor` function is smart and should not compute twice the `cor` function between col A and col B ; and between col B and col A ; thus 45 correlations should be computed)

To test the different genes (in a pairwise manner) :

cor(t(xx), method="pearson")

here a 58,000 x 58,000 matrix , = 3.364e+09 correlations (or 1,681,971,000 correlations if `cor` function is smart). That's why R crashes, it will take to long to compute so many correlations.

---

> Edit based on OP comments

Use the coefficent of variation : https://en.wikipedia.org/wiki/Coefficient_of_variation :

dat.coeff.var <- apply(dat,1,function(x){sd(x)/mean(x)})