cannin/biostars_qa huggingface dataset:
Question: Because of splicing, one variant may affect different transcripts and the effect can be different from one transcript to another. In the case where there are multiple transcripts associated to one gene, Annovar will output the information of all the transcripts.

However for the SIFT prediction (or any other SNP annotation tool - Polyphen2, Mutation Taster..), even if there are different transcripts involved by the variant, there's only one score/prediction.

For example:

C6orf25:NM_138273:exon3:c.C432G:p.F144L,C6orf25:NM_025260:exon4:c.C523G:p.R175G,C6orf25:NM_138272:exon4:c.C523G:p.R175G,C6orf25:NM_138277:exon4:c.C523G:p.R175G

Here there are 4 different transcripts listed, however the only SIFT prediction is `0.34/T`.

So my question is the following: how do I know which transcript was taken into consideration for the prediction? If there's only one prediction, does that mean that it applies to all the transcript?

I am trying to understand how I can interpret the SIFT prediction in that case and any help in that regard is welcome.

Answer: On the Annovar website, it's stated that

> ... it can have multiple scores such as "1.0;1.0;1.0;1.0;0.993;1.0;1.0;1.0;0.999" and multiple predictions such as "D;D;D;D;D;D;D;D;D", probably due to multiple transcriptional isoforms. In this case, only the largest score (1 is the largest score among multiple isoforms), as well as its associated D/P/B annotation, will be used in ANNOVAR database."

http://www.openbioinformatics.org/annovar/annovar_filter.html#ljb23

I believe the score that ANNOVAR reports is the most deleterious score out of all the transcripts. If you want scores that are transcript specific, or to figure out which transcript that score belongs to, I suggest either going back to dbNSFP, or trying out VEP or SnpEff.

Hope that helps
Question: Hi everyone!

I'm quite new to this genomics stuff, so sorry if this question seems silly (or just cheeky). I have looked for similar questions, but haven't found any (but I'm terrible at searching, though).

So. I have Illumina reads for a E.coli strain that has GFP inserted in it. I know roughly where the GFP should be, and I have the reference genome for the unmodified E.coli strain. I want to assemble a new reference for my modified strain.

From what I've seen so far, the idea would be to do a de novo assembly, then some scaffolding and mapping to the reference I have, using something like Mauve. But as my sequences map almost perfectly to the reference, and the only important breakpoint is the GFP insertion, can't I use this information to improve my assembly? If so, how?

Thanks!
Pablo

Answer: It's called reference guided assembly. Many assemblers can do it..
Question: I was wondering if there is any easy way to do an identity cutoff (e.g. Sequences down to 95% Id) from the blast output locally .I had a look at blastall usage options but haven't seen any flag (or parameter) to do it. There is e-value cutoff option but not %id or score. I have a perl script to do this cut off in the blast output but the script is huge almost 50-60 lines so trying to avoid that and want to use an easy script if any of you got any.

Answer: I think, I can do it this way: once I got the normal blast output then I can filter it.

awk '{OFS="\t"; if($%id_column_number>desired cutoff%)print$column_numbers_I_want_with_commas}' blast_output.file >> output

For example:

awk '{OFS="\t"; if($3>95%)print$1,$4,$5,$8}' blast_output.file >> output
Question: Hi,

I would like to use the "clump" function for my association analysis results which use 10000 Genomes imputed data (and I have many files), so I guess I need to use the 10000 Genomes genotypes as a reference file, but I am having trouble with this. I have downloaded the vcf files from the MACH website [here][1].

Then I thought I could use Plink 1.9 version that reads vcf and use directly the "clump" command here, but apparently there are issues with the vcf file for allele names (I would like to keep as many variants as possible), so I followed [this link][2] to convert from vcf to plink.

But I get a segfault error as soon as I type to this command:

bcftools annotate -Ob -x ID chr21.test.vcf

which says:

```
[W::vcf_parse] INFO 'NS' is not defined in the header, assuming Type=String
Encountered error, cannot proceed. Please check the error output above.
```

I have tried to use `bcftools reheader -h` but it doesn't work, can someone please help?

Thank you!!

[1]: http://www.sph.umich.edu/csg/abecasis/MACH/download/1000G.2012-03-14.html
[2]: http://apol1.blogspot.co.uk/2014/11/best-practice-for-converting-vcf-files.html

Answer: If you don't need to deal with split multi-allelic variants, Plink 1.9's `--set-missing-var-ids` flag should do a good enough job of assigning names. (You'll need a build from November 25th or later.)
Question: Dear all,

talking about the pathogenicity predictors on cancer mutations, what algorithms or meta-predictors would you recommend to use ? Among possible choices : CADD, MutationTaster, FATHMM, CHASM, Condel CanDrA , or any other predictors/meta-predictors.

thank you,

bogdan

Answer: Take your pick...

This is not a complete listing, as there are many more.

#Missense predictions

- SIFT
- Align-GVGD
- Mutation Assessor
- PANTHER
- MAPP
- MUpro
- MutationTasster
- I-Mutant
- KD4v
- FATHMM
- dbNSFP
- PolyPhen-2
- LS-SNP/PDB
- SNPeffect 4.0
- PMut
- SNAP
- PhD-SNP
- SNPs&GO
- Parepro
- CanPredict
- nsSNPAnalyzer
- MutPred2
- HANSA
- MutationTaster

#Splice predictions

- NNSplice
- GeneSplicer
- MaxEntScan
- HSF
- NetGene2

#Protein modelling (from amino acid sequence)

- Protein
Model Portal

*[uses various modelling algorithms and produces PDB files, which can be loaded into protein viewers like Jmol]*

#Non-coding (i.e. regulatory)

- CADD (germline
variants)
- DANN
(germline variants)
- FATHMM-MKL
(germline variants)
- GWAVA
(germline variants | somatic mutations)
- Funseq2 (somatic
mutations)
- SurfR
(rare variants | complex disease variants | all other variants)

# Other

- GWAS3D / GWAS4D

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

For further reading:

- UK perspective, refer to the guidelines by the Association for
Clinical Genomic Science:
http://www.acgs.uk.com/quality/best-practice-guidelines/
- US perspective, refer to American College of Medical Genetics:
https://www.acmg.net/docs/Standards_Guidelines_for_the_Interpretation_of_Sequence_Variants.pdf
Question: Hi everybody,
I'm trying to run this pipeline with bigwif files

computeMatrix reference-point --referencePoint TSS --beforeRegionStartLength 30 --afterRegionStartLength 300 -R transcribed_genes_sort.txt -S BJ_rep1_neg_strand.bw BJ_rep1_pos_strand.bw BJ_rep2_neg_strand.bw BJ_rep2_pos_strand.bw -bl DACblacklist.bed.gz --skipZeros -p 1 -o matrix_TSS.gz

and I get this error

[bwHdrRead] There was an error while reading in the header!
Traceback (most recent call last):
File "/home/projas/.local/bin/computeMatrix", line 11, in
main(args)
File "/home/projas/.local/lib/python3.5/site-packages/deeptools/computeMatrix.py", line 414, in main
hm.computeMatrix(scores_file_list, args.regionsFileName, parameters, blackListFileName=args.blackListFileName, verbose=args.verbose, allArgs=args)
File "/home/projas/.local/lib/python3.5/site-packages/deeptools/heatmapper.py", line 235, in computeMatrix
chromSizes, _ = getScorePerBigWigBin.getChromSizes(score_file_list)
File "/home/projas/.local/lib/python3.5/site-packages/deeptools/getScorePerBigWigBin.py", line 161, in getChromSizes
fh = pyBigWig.open(fname)
RuntimeError: Received an error during file opening!
However, the header of the bigwig files is something like:

chr1	13264	13314	1.8
chr1	13444	13454	3
chr1	13454	13464	4.6
chr1	13464	13474	8.6

Any ideas?
Thanks!

Answer: For those wondering, the cause of the problem was that the input files were actually bedGraph (text files) rather than bigWigs (binary format). Presumably there was something wrong with the conversion command with UCSC tools.
Question: Is there a way to generate a heatmap of a complex of transcription factors (in more detail, I am trying to generate a heatmap of the 7SK snRNP which is composed of Hexim, Larp7, and Cdk9) at the TSS. I can easily generate individual heatmaps for each of the separate TF's themselves using ngsplot or deepTools, but I was wondering if there was a way to combine all three into a single heatmap?

I also have some experience with R and Bioconductor so any packages can also be submitted as answers.

Thanks in advance

Answer: Just use the "release-1.6" branch of deepTools, you can then include multiple samples in a heatmap (or profile plot).
Question: I have RNA-seq data for 20 samples with 2 condition and 2 sex (male, female : control, treatment). I am very new to RNA-seq analysis and am trying to find the DEGs using DESeq2. Since I want to have the normalization to be calculated based on all samples, I get the rld based on all and then I will use contrast to find the DEG for 1. tratment vs. control male 2. tratment vs. control female.

For the QC, the PCA plot separates the male and female but the control and treatment is not separated very well on the PC2 for male C1 (Figure 1)![enter image description here][1]

Then I decided plot PCA only for male samples. again PC1 does not separate the control and treatment male samples but PC2 is kinda separating them. (Figure2)

![enter image description here][2]

My questions are:

1. What numbers on the PCA plot (x and y axis) decide the separation? I read in a [website][3] that samples that are at PC1>0 are outlier. is that true?

2. Can I just look at figure 1 and remove male c1 and continue DEG analysis with 9 samples? or I should definately plot figure 2 ?

3. If I need to consider Figure 2, can I rely on only the PC2 which separates control and treatment and continue the DEG analysis or I should remove C1, tr4, tr5 samples and then work on DEG analysis based on remaining 7 samples?

[1]: https://i.imgur.com/UgC0byK.png
[2]: https://i.imgur.com/mpnghro.png
[3]: https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#pca-and-sample-heatmaps

Answer: Hey,

Regarding the website to which you are referring, that was written by Wolfgang Huber. He is *not* saying, in his tutorial, that PC1>0 is a global definition of an outlier from a PCA bi-plot. It is just in *his* example that he has decided to select 2 outlier samples via a vertical cut-off line drawn at PC1>0. In another experiment, the cut-off for outliers may be at PC1 < -500, or PC1 > 42. The units on these axes are, well, essentially unitless.

Looking at your first plot, you cannot really remove male sample C1 just because it groups away from the main group of samples. You should first investigate *why* it may be grouping in the way that it is [grouping]. Is it a sample mix-up?; is it a male with Klinefelter Syndrome?; is the clinical information for this sample incorrect?

For your third question, you need to consider the fact that `sex` represents a substantial amount of variation in your dataset, i.e., it is greater than the very condition of interest that you are studying - was this expected?; did you include `sex` in your design formula as a covariate? If you include it as a covariate in your design formula, normalise your data again, and then run rld function with `blind=FALSE`, do you see the same PCA bi-plot?

You can also segregate your dataset based on `male` | `female` and then cross-compare results in a sort of meta-analysis, if you wish. `sex` is a clear and major source of variation among your samples. If you proceed with this route (sample segregation based on `sex`), then I see no reason to exclude c1.

Kevin
Question: Can EnhancedVolcano label genes that aren't significant or meet fold change cut off? Can we select? I tried to but it comes up as a geneid on the x-axis with vertical bar. Just wondering. 4.8 seems to suggest no.

How does EnhancedVolcano determine which genes to label? Is there a way to label more to the ones selected automatically? I guess following 4.8, you can selectLab but they have to be present already in lab. What happens if you selectLab genes that interfere with ones automatically selected by EnhancedVolcano, do they get bumped or removed?
Can we remove select labels we don't want? It seems no matter how genes are labeled, it picks the same ones to label.

The default is Log2FC=2, so means minimum 4 fold change? or total 2, so 2 fold change, log2 = ratio 1?

Should I use lfcShrink(res) to get volcano plot or is just res fine? The vignette seems to suggest lfcShrink but pdf manual example has just res. Not sure what the proper or standard way for volcano plots from RNA-seq data run through DESeq2 is and how people normally publish them.

This package does look great. I'm pretty beginner at this and this is pretty user-friendly with nice vignette.

Answer: I wrote *EnhancedVolcano* - thanks for your comments.

You mentioned section 4.8: 4.8 Only label key transcripts. This should also label genes that do not pass your thresholds. However, the points for these non-statistically significant genes is usually dense, at the bottom of the plot; so, labeling all of them in a clear way is almost always impossible. If you want to favour the labeling of certain genes, then you can (I believe) re-order your results data-frame to have these genes appear in the first few rows of your data-frame.

You can also supply custom labeling and colouring via 4.11 Override colouring scheme with custom key-value pairs:

For the fold-changes, you can technically supply any numbers as the value of `x`. These can be log2 fold changes, linear fold changes, etc. The corresponding cut-off, as `FCcutoff`, is not assuming any distribution. If you supply linear fold changes to `x` and have `FCcutoff=2`, then the cut-off is linear FC = 2. If you supply log2 fold-changes as x and have `FCcutoff=2`, then this cut-off relate to log2FC = 2

`lfcshrink()` can be used when you have initially used `betaPrior = FALSE` during the `DESeq()` function. `lfcshrink()` was a recent implementation into DESeq2, so, I am not sure of the exact situations where it should (or should not) be used.

Kevin
Question: I am trying to compile shapeit4 on my cluster without root permissions. The makefile is as follows:

#COMPILER MODE C++11
CXX=g++ -std=c++11

#HTSLIB LIBRARY [SPECIFY YOUR OWN PATHS]
HTSLIB_INC=/share/apps/genomics/htslib-1.10/include/
HTSLIB_LIB=/share/apps/genomics/htslib-1.10/lib/libhts.a

#BOOST IOSTREAM & PROGRAM_OPTION LIBRARIES [SPECIFY YOUR OWN PATHS]
BOOST_INC=/share/apps/boost-1.71.0/include/
BOOST_LIB_IO=/share/apps/boost-1.71.0/lib/libboost_iostreams.a
BOOST_LIB_PO=/share/apps/boost-1.71.0/lib/libboost_program_options.a

#HTSLIB LIBRARY [SPECIFY YOUR OWN PATHS]
#HTSLIB_INC=/software/UHTS/Analysis/samtools/1.4/include
#HTSLIB_LIB=/software/UHTS/Analysis/samtools/1.4/lib64/libhts.a

#BOOST IOSTREAM & PROGRAM_OPTION LIBRARIES [SPECIFY YOUR OWN PATHS]
#BOOST_INC=/software/include
#BOOST_LIB_IO=/software/lib64/libboost_iostreams.a
#BOOST_LIB_PO=/software/lib64/libboost_program_options.a

#COMPILER & LINKER FLAGS

#Best performance is achieved with this. Use it if running on the same plateform you're compiling, it's definitely worth it!
#CXXFLAG=-O3 -march=native
#Good performance and portable on most intel CPUs
CXXFLAG=-O3 -mavx2 -mfma
#Portable version without avx2 (much slower)
#CXXFLAG=-O3

LDFLAG=-O3

#DYNAMIC LIBRARIES
DYN_LIBS=-lz -lbz2 -lm -lpthread -llzma -lcurl -lcrypto

#SHAPEIT SOURCES & BINARY
BFILE=bin/shapeit4
HFILE=$(shell find src -name *.h)
CFILE=$(shell find src -name *.cpp)
OFILE=$(shell for file in `find src -name *.cpp`; do echo obj/$$(basename $$file .cpp).o; done)
VPATH=$(shell for file in `find src -name *.cpp`; do echo $$(dirname $$file); done)

#COMPILATION RULES
all: $(BFILE)

$(BFILE): $(OFILE)
$(CXX) $(LDFLAG) $^ $(HTSLIB_LIB) $(BOOST_LIB_IO) $(BOOST_LIB_PO) -o $@ $(DYN_LIBS)

obj/%.o: %.cpp $(HFILE)
$(CXX) $(CXXFLAG) -c $< -o $@ -Isrc -I$(HTSLIB_INC) -I$(BOOST_INC)

clean:
rm -f obj/*.o $(BFILE)

The compilation then fails on:

/usr/bin/ld: cannot find -lcurl

I understand that this means the compiler has looked in /usr/bin/ld and hasn't been able to find the lcurl library. Curl is installed on the cluster:

├── bin
│ ├── curl
│ └── curl-config
├── include
│ └── curl
├── lib
│ ├── libcurl.a
│ ├── libcurl.la
│ ├── libcurl.so -> libcurl.so.4.5.0
│ ├── libcurl.so.4 -> libcurl.so.4.5.0
│ ├── libcurl.so.4.5.0
│ └── pkgconfig
└── share
├── aclocal
└── man

I've tried to add several different lines to the makefile in order for it to find -lcurl, but I haven't succeeded so far. Can anyone advise me on what to add to the makefile?

Answer: Install it with BioConda https://bioconda.github.io/recipes/shapeit4/README.html
Question: Hi all,

I ran partitioned heritability analysis using LDSC and custom annotations and got negative heritability enrichment scores for 3 phenotypes out of 34. We were a bit confused about the negative heritability results and then I was advised to filter my summary statistics to keep only the SNPs which are present in the baseline model which I used to run partitioned heritability.

My first question is, how do you think keeping only the SNPs in my baseline would affect the heritability enrichments? I did a bit of googleing and found out that usually low heritability, sampling errors and low sample size drives negative heritability. How do you think filtering SNPs in my sumstats based on the baseline would change the heritability measure?

Secondly, a technical question: How can I filter my sumstats based on the baseline? Is there a function to use from one of the common tools or should I manually keep the SNPs that fall within genomic intervals in my baseline bed files?

Thanks!

Answer: For the record - I found out that my issue was too short annotations. I found out that your annotations must encompass at least 1% of the whole genome (1% of all SNPs in the 1KG phase 3 reference panel is around ~20k SNPs) in an old post at the LDSC google group.

However, based on my experience, you need at least around 90k-100k SNPs in your annotation for LDSC to work as intended and to minimize negative heritability estimates.

Hope this is helpful for others!
Question: Hi everyone,

I am trying to plot genomic regions using KaryoplotR. I have been following the tutorial page, I converted my data frame into Granges object, however the tool does not work with my data.
My data frame is simply a list of genomic location ['chr', 'start', 'end']

The following script doesn't give any warning or error message and only plots the chromosomes but not the genomic regions:

library(karyoploteR)
library(GenomicRanges)
library(data.table)

setwd("/path/to/dataframe")

df <- fread("for_ideogram.csv", header = T)

gr <- makeGRangesFromDataFrame(df, start.field="start", end.field= "end")

kp <- plotKaryotype(genome = "hg19", chromosomes="autosomal", plot.type = 1)

kpPlotRegions(kp, data=gr)

Does anyone know how to get this to work?

Thanks

Answer: Hi @gabi,

When karyoploteR creates the ideogram but nothing gets added to it with the other plotting functions, in many cases it is due to different chromosome names.

Are your gr chromosome names in UCSC style ("chr1", "chr2", ...) or in Ensembl/NCBI style ("1", "2", ...)? The default bioconductor chromosome names are UCSC style and these are the ones used by karyoploteR.

Please check if adding `seqlevelsStyle(gr) <- "UCSC"` before plotting solves the problem.

library(karyoploteR)
library(GenomicRanges)
library(data.table)

setwd("/path/to/dataframe")

df <- fread("for_ideogram.csv", header = T)

gr <- makeGRangesFromDataFrame(df, start.field="start", end.field= "end")
seqlevelsStyle(gr) <- "UCSC"

kp <- plotKaryotype(genome = "hg19", chromosomes="autosomal", plot.type = 1)
kpPlotRegions(kp, data=gr)

and I don't know the structure of your csv file, but if the first 3 columns contain the chr, the start and the end positions (like a bed file and with any names) you can use [`toGRanges`](https://rdrr.io/bioc/regioneR/man/toGRanges.html) simplify you code to this (toGRanges will change the chromosome names to the ones used by "hg19" in this case):

library(karyoploteR)

gr <- toGRanges("path/to/dataframe/for_ideogram.csv", genome="hg19")

kp <- plotKaryotype(genome = "hg19", chromosomes="autosomal", plot.type = 1)
kpPlotRegions(kp, data=gr)

If this does not work, I'd need to see the first couple of lines for you csv file to understand what's going on.
Question: Hi all,

I have some RNAseq data and I would like to take the results forward to build a network based on expression correlation (wgcna or similar).

I have several steps that I need to carry out:

- normalise data
- transform data
- subset data (not everything in the dataset is useful for the intended outcome so I need to extract only the useful samples)
- remove a known batch effect (I know I can model this for differential expression, but for building a network I think I need to remove it - correct me if I am wrong)

I am using DESeq2 for normalisation and transforming using the variance stabilising transform from the same package (as recommend in the wgcna manual).

I have noticed that the outputs of my exploratory analyses change depending on the order in which I carry out these steps, particularly PCAs. In most cases, the gross patterns in the data remain intact, but in some cases this is not true. My question is, what is the correct order in which to carry out these steps and why?

Currently, I am loading all of the available samples, normalising, transforming the normalised counts, removing the batch effect from the transformed data, then extracting the samples of interest. General discussion about the order in which these processes should be carried out is welcome, but specifically:

1. Would it be more sensible to extract the samples I am interested in first, and then run the downstream steps only on the samples I am interested in? I imagine this would affect the output as the geometric mean across the samples would change.
2. I am currently removing the batch effect from the vst-transformed data. Would I be better to remove the batch effect first from the normalised counts and then transform the batch-corrected data?

Thanks for your help
Question: Hi all,

I have **several BAM files which I would like to work in R environment**. This are big files with **more than 100M reads** each. Is there any way to convert them in GRanges object and then save it as RData compressed file?

So far what I've tried is the following R code:

library(GenomicRanges)
library(Rsamtools)

args = commandArgs(TRUE)
filename <- args[1]
name <- args[2]

## Create GRanges object from the BAM file and save it to disk
param <- ScanBamParam(what=c("qname","flag"))
b <- readGAlignments(filename, format="BAM", param=param);
save(b,file=paste(name,".RData",sep=""), compression_level=9);

This code gave me error messages because R is not able to allocate vectors more than 800 ~ 1000 Mb size, which I understand, but is there **any way to create GRanges object of large BAM files**?

Thanks for your help!!

Answer: It works, you just need a lot of memory and coffee. I have tested with a 21GB BAM file containing ~400Million aligned reads. The process grows to ~70GB in size and the loading takes several hours. Even though it works in principle, it might be a good idea to a basic summarization outside of R, e.g. read counting or coverage calculation. To save a little space you could maybe ommit the read names. Another possibility might be to split the input by chromosome, but then again, maybe R is not the the optimal software for this volume.