cannin/biostars_qa huggingface dataset:
Question: I have been working with TCGA cancer data to examine expression (RNAseqV2) and methylation (Illumina 450k) data. I want to look at sequencing data, but I'm a bit lost to what sort of information is available through TCGA. I want to examine whether there are nonsense mutations between positions 2000-2500 across all cancer types available in TCGA. What sort of resources/workflow should I expect?

Answer: If you are asking whether you can search for non-sense mutations between amino acids 2000-2500 of a particular gene in the TCGA sets, I would suggest using the Mutation Annotation Files (MAFs). These files have already been run through somatic variation calling so you won't have to deal with the sequencing data directly. You can filter out your gene of interest across all the MAF's in perl,awk,grep etc... Then in the `amino_acid_change` column you can search for integer values between p.2000-p.2500, and `trv_type` "Non-sense mutation".

Hope this helps.

Nick
Question: I'm trying to get a regex to work with rename; I've tried the approach of similar answered questions here but couldn't get the results I wanted.

The files are named as such:

SR1_S90_L001_R1_001.fastq.gz
SR1_S90_L001_R2_001.fastq.gz
Rinc_S96_L001_R1_001.fastq.gz
Rinc_S96_L001_R2_001.fastq.gz

And I would like to retain only the information prior to the first underscore and the _R1_ or _R2_ tags, like this:

SR1_R1_.fastq.gz
SR1_R2_.fastq.gz
Rinc_R1_.fastq.gz
Rinc_R2_.fastq.gz

Thanks in advance!

Answer: Try safe-batch-rename tool `brename` ( https://github.com/shenwei356/brename )

brename -p '^(\w+?)_.+_(R[12])_.+' -r '${1}_$2.fq.gz' # updated

# original answer
# brename -p '^(\w+)_.+_(R[12])_.+' -r '${1}_$2.fq.gz'
# if you have ran this, you can run 'brename -u' to undo.
Question: I want to determine the substitution error rate and the indel error rate for a given BAM file.

I've been reading over the following article.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/#B1

First off I just want to confirm my definitions of a substitution error and indel error are correct.

A substitution error would be when the sequencer substitutes a different base than the actual base in the sequence being sequenced. So if the actual sequence was ..ATGG.. but the sequencer read ..ACGG.., then the reading of C instead of a T would be a substitution error, since one of the reads has a C which shouldn't be there.

An indel error would be when the sequencer deletes or inserts a base that is different from the actual sequence being looked at. So if the actual sequence was ..GATG.. and the sequencer read ..GTG.. then the A being deleted would be an indel error., since the read is missing an A that should be there.

Assuming those are right, I am also confused with how to determine them from a BAM file. If a tool exists then that would be great

For substitution I know it would involve the quality score for a given base and the number of mismatches. By mismatches I just mean if a given base has a coverage of 100X and 98 are T but the other two aren't, then there are 2 mismatches for that position. I'm just not sure how exactly I would combine these to find a rate.

For indel errors I know it would involve homopolymers like mentioned in the paper but I have no idea how to find the rate.

Answer: BBMap gives very complete statistics for substitution, insertion, deletion, and N (no-call) rates; they get printed to the screen after mapping. Also, it can generate a histogram with the `mhist` flag, showing the match/sub/ins/del/N rate by position in the read; `ehist` shows how many reads have a specific number of substitutions; and `indelhist` gives the counts of indels of specific lengths.

bbmap.sh ref=reference.fa in=reads.fq out=mapped.sam mhist=mhist.txt ehist=ehist.txt indelhist=indelhist.txt

These histograms can also be generated by Reformat from a sam or bam file, but only if the cigar strings are in sam 1.4+ format (with X and = in the cigar strings instead of M). Most aligners do not generate that kind of cigar string, but you can easily check by looking at one of the mapped reads in the bam file.

If you want to exclude actual polymorphisms in the genome that do not come from sequencing, it's probably best to call variations and mutate the reference according to the calls, then map again. Alternately, for SNPs at least, you can get an estimate of the number of substitutions by running an error-correction program which will generally tell you how many errors it finds and corrects.
Question: I checked inbreeding coefficient (F) on my samples (around 200) using plink

plink --file mydata --het

and found the distribution of F pretty symmetric:

and the short summary:

> summary(het$F)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.35990 -0.07199 0.01568 0.02755 0.13460 0.40790

The F values are quite big (either positive or negative), which made me worry about sample quality. Presumably they are all unrelated. Can anyone tell me whether this distribution of F indicates some problems and what they are? Or, is there a rule of thumb about what F value is considered normal for unrelated, cleaned samples? Thanks!

Answer: For relatedness try IBS/IBD estimation:

- plink 1.9 https://www.cog-genomics.org/plink/1.9/ibd
- old website: http://zzz.bwh.harvard.edu/plink/ibdibs.shtml

PI_HAT:

- Identical twins, and duplicates, are 100%identical by descent (Pihat 1.0)
- First-degree relatives are 50% IBD (Pihat 0.5)
- Second-degree relatives are 25% IBD (Pihat 0.25)
- Third-degree relatives are 12.5% equal IBD (Pihat 0.125).
Question: Hi.

I am wondering a term with SIFT score.

I think that SIFT refers to some measurement of SNPs, and while reading Annovar paper, I saw below sentence as follows:

> Finally, Annovar can filter specific variants such as SNPs with >1% frequency in the 1000 Genomes Projects, or non-synonymous SNPs with SIFT scores > 0.05.

Regarding above sentence, I ask you two questions.

1. I think that 1% frequency is a little bit low allele frequency. Dose it have an effect to filtering irrelevant snp variants? I don't think so..
2. SIFT-score threshold is about 0.05 as shown in above sentence. What does SIFT means about and threshold of 0.05 might be effect on filtering variants?

Answer: To answer your first question, 1% is the standard cutoff used to describe the difference between "common" and "rare" variants. Depending on your study, you might want to change that. For example, in a GWAS for a common trait, you might be interested only in variants that are above a certain frequency in the population, whereas if you're looking at rare Mendelian traits you might only want very low frequency variants. You may also want to narrow this down to a specific population, eg for a GWAS in African Americans, you would be interested in variants common in African populations. Steve mentioned the VEP, which allows you to filter variants by frequency, choosing your own frequency, > or < and pick a population of interest.
Question: Hi friends

I installed Tophat2 with all dependencies (bowtie2, Boost, Samtools) but when I want to run tophat, it shows this error:

```
File "/home/tophat-2.0.14/bin/tophat", line 1003
except getopt.error, msg: ^
SyntaxError: invalid syntax
```

My python version is 3.4 , how can I fix this problem?

Thanks

Answer: Try installing python2.
Question: Trying to use [tbl2asn](http://www.ncbi.nlm.nih.gov/genbank/tbl2asn2/)

The directions are [here](http://www.ncbi.nlm.nih.gov/books/NBK53709/#gbankquickstart.Submission_using_tbl2asn)

It starts on Step C where you download the tbl2asn file and rename it. It says "set permissions as required for your operating system." but i have no idea how to do that.

I moved it into a new folder in my applications and then:

export PATH=/Applications/Command_line_apps:$PATH

in the terminal . tried to call it using tbl2asn in the command prompt but I got back

bash: /Applications/Command_line_apps/tbl2asn: Permission denied

How do fix the permissions?

Answer: Hello jolespin - this post does not fit the main topic of this site, as is not a bioinformatics question. I'd suggest that you do some basic reading on file permissions in UNIX.

FWIW, 'chmod +x yourfile' should make it executable, which may solve your problem.
Question: Hi,
Is there a tool out there which will allow me, in any way, to determine how similar two fastq file sets (paired-end) are? It could be any metric like number of identical reads etc. or any other metric which can be relevant in this case.

I need this to diagnose the reason behind low agreement of variant calls between two identical runs: if the fastqs are quite similar to each other, then it was the variant-calling pipeline and not the upstream bench-work.

Thanks!

Answer: Maybe you could try [commet][1]? It was designed for metagenomics, but it allows you to compute a distance between two fastq files.

[1]: https://colibread.inria.fr/software/commet/
Question: Dear All,

I would like to know about few docking softwares, freely available ones.

Thanks

Answer: Honestly, depends on which *docking* you are talking about. The mechanics of the software (i.e. the algorithms, scoring functions) depend on its purpose. If you want to do protein-protein docking, Autodock might not be the best option. For protein-ligand docking, it might. You also have protein-nucleic acid, protein-glycan, multiple proteins, etc.

Autodock is a good option for protein-ligand docking. If you want an overview for protein-protein/DNA/etc docking, have a look at these reviews: [1][1], [2][2]. Nearly all that software is free for academic use.

[1]: http://www.ncbi.nlm.nih.gov/pubmed/19462412
[2]: http://onlinelibrary.wiley.com/doi/10.1111/febs.12771/abstract
Question: Hi!

I have assembled my PacBio FASTQ reads with Canu. Now I would like to polish\correct the assembly by mapping these PacBio FASTQ reads on the assembly itself. I heard about Quiver in the SMRTanalysis pack, but I'm wondering if there is any alternative software.

Thanking you in advance for your help!

Answer: The tools [Pilon][1] and [Icorn2][2] performs assembly polishing with Illumina data. I am not sure if Icorn2 also takes PacBio reads for correction. Running additional rounds of Quiver may be the helpful to use the PacBio reads efficiently.

[1]: http://software.broadinstitute.org/software/pilon/
[2]: http://icorn.sourceforge.net/
Question: I have illumina paired-end whole-genome sequencing reads which I have map to around -400 reference plastid genomes. After getting mapped reads, I have to assemble as de novo plastid genome.

1. Do I have map reads to invidiual reference genomes one by one, or can I download all genomes at one go and index as one reference genome?. Do bwa or bowtie has

enough memory to index 400 genomes as one reference index genome?

2. Do you think which one is best method?. Mapping individual genome or all genomes indexed as one?

3. If I have to map individually, can I combine all bam file together and Can I convert to fastq file using bam2fastq tool (in picard) for denovo assembly?

Answer: I don't fully understand your introductory sentence as stated, but I can answer your questions to the best of my understanding:

1. You can concatenate all of your genomes into one large file and then index that composite genome. Just make sure your naming convention for each component in the FASTA is logical so that you can understand your results downstream. bwa and bowtie don't "have their own memory," but if you're using a 64-bit system you shouldn't run into any size limitation troubles, especially with plastid sequences.

2. **EDIT:** I say map to a single composite reference. You can pass parameters into bowtie to limit your mappings on the front end, thus saving computational time and making it easier to isolate the most accurate mappings for each read.

3. You can merge BAM files using samtools merge, but of course you wouldn't need to if you proceed according to my recommendation. There are a number of tools from converting from BAM back to FASTQ, and the Picard tool should work just fine. It does have trouble with paired-end mappings in certain circumstances though, and if you run into trouble using Picard I'd suggest bedtools bamtofastq.

Is the idea here that you're going to map to a bunch of plastid reference sequences from various organisms, and then convert the aggregate mappings back to FASTQ and perform an assembly from them? If so, I say concatenate the reference sequences into one file and map against that. That way it will be easier to retain the best mappings up front, especially if you don't care about *which* reference you're mapping to.
Question: Hello all,

`bcftools` have these logical operators that can be used in filtering expressions:

&& (same as &), ||, |

What's the difference between `||` and `|`? Can someone provide an example and/or usecase for clarification?

The [manual][1] have this example:

QUAL>10 | FMT/GQ>10 .. true for sites with QUAL>10 or a sample with GQ>10, but selects only samples with GQ>10
QUAL>10 || FMT/GQ>10 .. true for sites with QUAL>10 or a sample with GQ>10, plus selects all samples at such sites

But this is not clear to me.

fin swimmer

[1]: http://www.htslib.org/doc/bcftools.html#expressions

Answer: Say your VCF contains the per-sample depth and genotype quality annotations and you want to include only sites where one or more samples have big enough coverage (`DP>10`) and genotype quality (`GQ>20`). The expression `-i 'FMT/DP>10 & FMT/GQ>20'` selects sites where the conditions are satisfied within the same sample:

bcftools query -i'FMT/DP>10 & FMT/GQ>20' -f'%POS[\t%SAMPLE:DP=%DP GQ=%GQ]\n' file.bcf

49979 SampleA:DP=10 GQ=50 SampleB:DP=20 GQ=40

On the other hand, if you need to include sites where both conditions met but not necessarily in the same sample, use the && operator rather than &:

bcftools query -i'FMT/DP>10 && FMT/GQ>20' -f'%POS[\t%SAMPLE:DP=%DP GQ=%GQ]\n' file.bcf

31771 SampleA:DP=10 GQ=50 SampleB:DP=40 GQ=20
49979 SampleA:DP=10 GQ=50 SampleB:DP=20 GQ=40

This example is taken from http://samtools.github.io/bcftools/howtos/filtering.html

---
**EDIT: (inserted by a mod)**

*Answer given on [github][1]:*

Well, sorry to demonstrate the difference on `&` and `&&` instead of `|` of `||`, but it's the same priniciple.

The manual page says it all:

QUAL>10 | FMT/GQ>10 .. true for sites with QUAL>10 or a sample with GQ>10, but selects only samples with GQ>10
QUAL>10 || FMT/GQ>10 .. true for sites with QUAL>10 or a sample with GQ>10, plus selects all samples at such sites

Or you can try to run yourself:

$ bcftools query -f'[%POS %SAMPLE %DP\n]\n' -i'FMT/DP=19 | FMT/DP="."' test/view.filter.vcf
3162006 A 19

3162007 A .
3162007 B .

$ bcftools query -f'[%POS %SAMPLE %DP\n]\n' -i'FMT/DP=19 || FMT/DP="."' test/view.filter.vcf
3162006 A 19
3162006 B 1

3162007 A .
3162007 B .

[1]: https://github.com/samtools/bcftools/issues/856#issuecomment-416140784
Question: I've downloaded a fastq file from SRA (http://trace.ncbi.nlm.nih.gov/Traces/sra/) containing reads from a paired-end Illumina 101 bp RNAseq experiment. The only problem is, it contains both read pairs in a single file, whereas I need separate files with all the _1.fq reads in one and the _2.fq reads in another.

Can anybody help? I'm aware of the fastq-dump tool within the SRA Toolkit, but I couldn't get it to work when I was originally downloading the data.

Many thanks in advance.

My fastq file looks like this:

$ head sra_data.fastq

@SRR1659960.1.1 1 length=101

NAGAAATGAATGAGCCTACAGATGATAGGATGTTTCATGTGGTGTATGCATCGGGGTAGTCCGAGTAACGTCGGGGCATTCCGGATAGGCCGAGAAAGTGT

+SRR1659960.1.1 1 length=101

#1=BDDDDDHFBFIEHHHHAG<HE@HGGE@HHFGHGGHHFHIHG@FFGGGHIIIIIFAC=F@GEGEECCDCECCBBBBCCCD>9599>C:@>5@9>?CCCD

@SRR1659960.1.2 1 length=101

CCCACTTCCACTATGTCCTATCAATAGGAGCTGTATTTGCCATCATAGGAGGCTTCATTCACTGATTTCCCCTATTCTCAGGCTACACCCTAGACCAAACC

+SRR1659960.1.2 1 length=101

<7?BD?DD<DFFABBEHEEFHII>C:BCDD?<C?FFC4E>@DEF>?FGHDFBBCG8??DGGIII:BF@C=FFC;C=D;@?EA76?DDBEC?>>ACCCABBB

@SRR1659960.2.1 2 length=101

NATAAAGTGTATGACAAATATACAAGGCTCCTAATATTGGTTTAACTTGGAGAAGTAGGTAAAGGAAGAAGGGNAAAGGAAATAGACAAAAAGACTACAGT

Answer: Use Reformat from the BBMap package:

reformat.sh in=sra_data.fastq out1=r1.fq out2=r2.fq interleaved
Question: Dear all,

has anyone done any benchmarking on speeding up long read alignment algorithms ?

I mainly use minimap2, but its' runtime varies by a factor of 10 across our cluster. I've been trying mm2-fast https://github.com/bwa-mem2/mm2-fast, the partially accelerated version, but without much success so far.

Is for example PAF output faster than SAM ?

Have others worked out how to scale minimap2 for Promethion scale datasets ? I expect LRA https://github.com/ChaissonLab/LRA has a similar runtime from their presented results, and others seem slower still (ngmlr etc).

Thanks

Answer: 1) Yes, when Minimap2 makes paf files it is faster than when it makes sam files.

2) For a speedup at cost of accuracy you can increase the minimizer length ("-k") and window length ("-w").

3) You can increase "-I". If the reference is larger than 4 Gbp, this will accelerate Minimap2 at cost of increased RAM consumption.

Also, see https://github.com/lh3/minimap2/issues/322
Question: I am interested in assessing the clonality of a tumor sample I have which is also paired with a normal tissue control from the same patient (~30x coverage). Specifically, I would like to know if there are any significant sub-clones present within the tumor sample.

I have come across numerous tools which look like they may be able to help me, including: [cloneHD][1], [SubcloneSeeker][2], [ABSOLUTE][3], [CITUP][4], [THetA2][5], [PhyloWGS][6] .. others?

Does anybody have an experience using any of these tools, or any thoughts on which one I should try out?

[1]: http://www.sciencedirect.com/science/article/pii/S2211124714003738
[2]: http://www.genomebiology.com/2014/15/8/443
[3]: http://www.nature.com/nbt/journal/v30/n5/full/nbt.2203.html
[4]: http://bioinformatics.oxfordjournals.org/content/31/9/1349.full
[5]: http://bioinformatics.oxfordjournals.org/content/30/24/3532.full#ref-6
[6]: http://www.genomebiology.com/2015/16/1/35

Answer: We've been working on integrating clonality and heterogeneity estimation tools into bcbio (https://github.com/chapmanb/bcbio-nextgen). It's still a work in progress but we've been evaluating against internal datasets where we have external predictions of normal contamination. We've had the most success with:

- Battenberg from Sanger (https://github.com/cancerit/cgpBattenberg). This calls CNVs and provides estimates of normal contamination. It's also a required input for PhyloWGS.
- PhyloWGS. We don't currently have a way to validate heterogeneity but get a useful set of trees to estimate how much of a mix is in a tumor sample.
- BubbleTree from MedImmune (http://www.bioconductor.org/packages/release/bioc/html/BubbleTree.html) which provides estimates of normal contamination and also clonality. We haven't dug much into the heterogeneity estimates of this yet.

We also looked at THetA2, but didn't have good luck with it, so spent more time with the above three tools.

Hope this helps provide some useful directions.
Question: Hello,

I am trying to use bedtools subtract in the following way:

$ bedtools subtract -a A.bed -b B.bed

but I get the following error:

ERROR: file A.bed has non positional records, which are only valid for the groupBy tool

When I tried to Google the answer, none of the solutions provided seemed to help. Here is the head of A.bed:

$ head A.bed
chr1 1 249250621
chr2 1 243199373
chr3 1 198022430
chr4 1 191154276
chr5 1 180915260
chr6 1 171115067
chr7 1 159138663
chrX 1 155270560
chr8 1 146364022
chr9 1 141213431

The cross checked that the columns are tab delimited and the three important columns are in the right order. Additionally, there are no invisible characters after the final column that might throw off Bedtools.

Any insight into the problem would be greatly appreciated!

Answer: I found the solution. I made A.bed in Windows and B.bed in Linux, so the two files had conflicting end characters. When I turned the Windows end characters into Linux end characters, Bedtools started working properly again
Question: I have data set between two comparison or two condition such as HSC and CMP , so i want to do a scatter plot with along with the regression value. Meanwhile when i m trying to label the sample they are all labelled into same color , i want different color for the labelling so that it can be distinguished

Here is my sample data set

gene HSC CMP
ENSG00000158292.6 1.8102636 2.456869
ENSG00000162496.6 2.6796705 6.203838
ENSG00000117115.10 3.4509115 5.555739
ENSG00000159423.14 3.6809277 5.063446
ENSG00000053372.4 5.7089974 6.851090
ENSG00000127423.8 4.4894292 5.996304
ENSG00000242125.3 10.6258802 11.715932

Here is my code

library("ggpubr")
ggscatter(data1, x = "HSC", y = "CMP",
add = "reg.line", conf.int = FALSE,
cor.coef = TRUE, cor.method = "pearson",
color = "black", size=1)

Any suggestion or help would be highly appreciated

Answer: If I understand well you want a scatterplot with a linear regression ?

>Edit to take your comment into account i.e. the color variable

using ggplot if data1 is your dataframe :

ggplot(data1,aes(x=HSC,y=CMP,col=variable)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE)
Question: Hello

I have WIG file and need to find genes regulated by the protein which is alternatively spliced (to be precise - the type of alternative splicing should be a5ss - alternative donor sites)

Please advise what programs should I use for this purpose.

Answer: [MATS][1] is my favorite tool for splicing event analysis

[1]: http://rnaseq-mats.sourceforge.net/

As mentioned in the previous comment, MISO is another option. In both cases, you need `.bam` files for RNA-Seq data (neither uses ChIP-Seq data, which is probably not a reliable tool for analyzing splicing events).
Question: Hello.
I am trying to run a blast against NR locally, and I am looking to find a way to run it as fast as I can. I recalled from my memory that I used LC_ALL=C fgrep -wf- for a "quicker" search file from one file to another, so I'm wondering if I can do something similar with the blast.

In the same line with fgrep I tried

blastp -query FILE -db LC_ALL=C /home/Protein_DataBases/nr/nr.fasta -evalue 1e-5 -outfmt 6 -num_threads 25 but it doesn't work. Does anyone knows if I can do this action??

Thank you in advance!

Answer: Don't put `LC_ALL=C` as a parameter to blast! That won't work.

You can set LC_ALL=C before a command:

LC_ALL=FOO; echo $LC_ALL

prints:

FOO

it gets applied to all commands run in that shell.

The best recommendation would be to export the variable in bash (preferable upon initialization) so that it is always applied.

export LC_ALL=C

you may run into various problems otherwise when it comes to bioinformatics tools and processes. Specifically sorting will be byte-wise when is set to `C` and alphabetical otherwise.
Question: How can I best plot a histogram for billions of genotype quality values?

I have a simple one column file with billions of genotype quality values. The file is several GB uncompressed.

Is there a statistics library in Python or R that can build up a histogram by streaming trough the data?
Instead of loading everything to memory and then creating the histogram? I prefer using all of the data versus sampling it.

Or do I have to write a script first to collect the counts per bin and then give those count per bin to R for plotting? This functionality feels like it should already exist in a stats library some where.

I know the min and max of the values and would be able to specify a bin size.

Answer: Meet the selling point of Datashader: http://datashader.readthedocs.org/en/latest (even though I tend to filter and subset the data whenever possible).