cannin/biostars_qa huggingface dataset:
Question: I have some problems with understanding of the concepts of SNP and related MAF, Alleles and NGVS Names.

I was looking at rs2476601. From "Allele" Information I could understand that

1. Reference Allele is A and variation results in G; and
2. A is the minor allele and has a frequency of 2.7%

Does it mean that A, our reference allele, is a minor allele and G, our variation, is a major allele? Why is then A called a 'reference allele' if it is a minor allele? Wouldn't it be logical to call a major allele, in our case G, to be a reference allele? Or a reference allele is the allele in the genome version used in mapping etc?

From 'HGVS Names':

- NC_000001.10:g.114377568A>G
- NC_000001.11:g.113834946A>G
- NG_011432.1:g.41808C=
- NG_011432.1:g.41808C>T
- NM_001193431.1:c.1858C=
- NM_001193431.1:c.1858C>T

There are sequences where A is mutated to G, but also sequences where C is mutated to T. I do not understand how it is possible. If we are looking at sense strand at see an 'A' there, which changed somewhere to 'G', it means that on antisense strand we will see a 'T' which changed to 'C'. So, it must be T>C. Why they show C>T then?

Thank you in advance.

Answer: The reference allele is always the allele that is in the reference genome. The reference genome is made up of contigs of sequence, all of which were taken from individuals, so the reference at any given position is that allele that that person happened to have at that locus - at various loci that allele happens to be the minor allele. Indeed, for any given variant, the probability that the reference contains the minor allele is the MAF. Long term, the GRC, who produce the reference genome, are attempting to switch the allele in the reference to whatever the major allele is at that locus, but this is very much in the long term.

G/A vs C/T is down to people reporting the alleles on the opposite strand. Usually this occurs when someone is reporting against a wild-type gene, which will often have the major allele. They will report the variant as major/minor on the reverse strand and submit this to the database, however, as we saw, the ref/alt might actually be minor/major.
Question: Hi pals, I got some problems here.

I have an excel file with several genomic coordinates in the format "strand:start:end", and I also have access to my genome of interest in different formats (fna, gbk...). Now, what I want to do is to extract and obtain a file with all the sequences corresponding to my coordinates, I don't mind if the sequences have a header or a custom ID, I just want the sequences in the same order as I have them in my excel file. Is there any tool that can do this? It would be much appreciated since I have to extract more than 2000 sequences, and doing this manually is humanly impossible.

The next step is to perform a BLAST with another genome, also accessible in several formats. I'm using the blastn NCBI program and this script `-task blastn -query "query" -subject "subject" -oytfmt 6 >name`, how can I make the output file to only contain the most relevant matches? (E.g. only the matches that have an E value lower than 0.0001)

Thanks beforehand!

Answer: Save your excel data as tab delimited file.

Use [bedtools getfasta][1] to fetch sequences:

bedtools getfasta [OPTIONS] -fi -bed -fo

[1]: https://bedtools.readthedocs.org/en/latest/content/tools/getfasta.html
Question: It is quite common to find tools that report on statistics per variant/marker, given a vcf file, but does anyone know of a tool that can report on stats per sample?

More specifically, I am looking for a tool that can tell me the proportion of variants/markers that are heterozygous (1/0 or 0/1) for each SAMPLE, given a vcf file.

Answer: [RTG Tools][1] includes a `vcfstats` command that outputs basic statistics for every sample (or you can request specific samples). e.g. on a small simulated vcf:

```
$ rtg vcfstats family.vcf.gz

Location : family.vcf.gz
Failed Filters : 0
Passed Filters : 144

Sample Name: sm_mom
SNPs : 91
MNPs : 1
Insertions : 5
Deletions : 2
Indels : 0
Same as reference : 1
Missing Genotype : 44
SNP Transitions/Transversions: 1.74 (73/42)
Total Het/Hom ratio : 2.96 (74/25)
SNP Het/Hom ratio : 2.79 (67/24)
MNP Het/Hom ratio : - (1/0)
Insertion Het/Hom ratio : 4.00 (4/1)
Deletion Het/Hom ratio : - (2/0)
Indel Het/Hom ratio : - (0/0)
Insertion/Deletion ratio : 2.50 (5/2)
Indel/SNP+MNP ratio : 0.08 (7/92)

Sample Name: sm_dad
SNPs : 73
MNPs : 2
Insertions : 2
Deletions : 3
Indels : 0
Same as reference : 1
Missing Genotype : 63
SNP Transitions/Transversions: 1.87 (58/31)
Total Haploid : 19
Haploid SNPs : 17
Haploid MNPs : 0
Haploid Insertions : 1
Haploid Deletions : 1
Haploid Indels : 0
Total Het/Hom ratio : 2.59 (44/17)
SNP Het/Hom ratio : 2.50 (40/16)
MNP Het/Hom ratio : 1.00 (1/1)
Insertion Het/Hom ratio : - (1/0)
Deletion Het/Hom ratio : - (2/0)
Indel Het/Hom ratio : - (0/0)
Insertion/Deletion ratio : 0.67 (2/3)
Indel/SNP+MNP ratio : 0.07 (5/75)

[...]
```

[1]: http://realtimegenomics.com/products/rtg-tools/
Question: I recently encountered the term "BAF" in a post on Biostars. I've looked up the definition, which is "B-Allele Frequency" (the extended definition is no more enlightening; one allele is A and the other is B). To me, that has no meaning other than being "Allele Frequency" with perhaps implicit assumptions about ploidy that are not generalizable, or maybe assumptions about the major/minor allele counts (from sequencing depth), or perhaps population frequencies, or maybe the reference allele. Can anyone explain where BAF is a useful term, or when one should use "BAF" instead of "AF"?

Answer: I've mostly seen BAF used in the context of SNP arrays, when there are two probes (an A probe and a B probe) covering a specific position. I think the A probe is generally the reference sequence (or was at the time the array was designed), so BAF and AF end up being the same, with the exception that SNP arrays are obviously limited in what they can detect.
Question: Hi all,

We're performing NGS on cancer DNA samples with a custom panel of target regions in different genes.

From your experience & knowledge, which tools do you recommend using for our purpose:

1. Finding structural variants (SVs) such as inversions, large indels, translocations, tandem duplications, etc.
2. CNV (copy number variantion) discovery for only a **single** sample (many tools I've found were for multiple samples).

If I'm not mistaken I should look for tools that are suitable for exome data rather than whole genome because of our sequencing methods and custom panel.

Thank you very much in advance!

---

I realize this has been somewhat asked before but the posts I've found were at least 2-3 years old and didn't elaborate on cancer genomics or exome sequencing. It was also unclear to me which CNV discovery tools are suitable for single samples. I hope this isn't out of line.

Answer: Lumpy, Delly, Manta, and Pindel are good for detecting structural variation and will work on exome/targeted samples. For CNV detection, some recent review articles give an overview of what's available for your application:

http://bib.oxfordjournals.org/cgi/doi/10.1093/bib/bbu004

http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-732

These articles are not comprehensive. A couple of tools not mentioned there but worth checking out are CopywriteR (Bioconductor) and CNVkit (mine; Unix command line) -- these are relevant if you're using hybrid capture, but less so for targeted amplicon capture. In the latter case, try OncoCNV.

In general you can get a good list of tools that are worth using by looking at existing validated pipelines like bcbio-nextgen and speedseq. Even if you don't use these pipelines directly, you can learn from others' experience by reading the documentation they provide.
Question: Hello!
I'm trying to understand what is the best algorithm for GWAS nowadays.
I know we have many tools available like Plink and Hail, but currently, what is the best algorithm if I won't use any them? Let's say, write down a script in R or Python from scratch. Which statistical algorithm should I use?
Is it linear mixed models (LMMs)? I'm confused as we can have binary phenotypes (case/control) or quantitative phenotypes. LMM seems to address quantitative ones, but can it be used for case/control as well?
Actually, what is the state of the art for both/each of them? Pair-reviewed papers as references will be appreciated.
Thanks!

Answer: The main regression executed by Plink was introduced by EIGENSTRAT in ~2006; see https://www.nature.com/articles/ng1847 . This is actually straightforward to write in R/Python from scratch; the harder part is optimizing the implementation for large datasets.

The Firth regression added to Plink 2.0 to improve handling of rare variants and imbalanced binary phenotypes was motivated by https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4049324/ .

Mixed linear models provide better statistical power when you have lots of close relatives in your dataset, but are much trickier to solve; actually, this is still a significant research area. Two tools covering parts of the current state-of-the-art are SAIGE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6119127/ ; handles imbalanced binary phenotypes, but relatively slow) and fastGWA (https://www.nature.com/articles/s41588-019-0530-8 ; great speed, but doesn't support dosage data yet and uses a misspecified model for binary phenotypes).
Question: Hello Biostars!

I hope you are shining bright!

Today, I'm dealing with a very basic problem that's plotting lines color change. It's very easy in general plot function and ggplot.
But in case of a function by ChIPseeker - 'plotAvgProf2', I'm not able to figure how to give arguments as there's not default argument like 'col' or 'colour' in the function :

plotAvgProf2 package:ChIPseeker R Documentation

plotAvgProf

Description:

plot the profile of peaks that align to flank sequences of TSS

Usage:

plotAvgProf2(peak, weightCol = NULL, TxDb = NULL, upstream = 1000,
downstream = 1000, xlab = "Genomic Region (5'->3')",
ylab = "Read Count Frequency", conf, facet = "none", free_y = TRUE,
verbose = TRUE, ...)

Arguments:

peak: peak file or GRanges object

weightCol: column name of weight

TxDb: TxDb object

upstream: upstream position

downstream: downstream position

xlab: xlab

ylab: ylab

conf: confidence interval

facet: one of 'none', 'row' and 'column'

free_y: if TRUE, y will be scaled by AvgProf

verbose: print message or not

...: additional parameter

I don't even how to edit the function too :

> plotAvgProf2
function (peak, weightCol = NULL, TxDb = NULL, upstream = 1000,
downstream = 1000, xlab = "Genomic Region (5'->3')", ylab = "Read Count Frequency",
conf, facet = "none", free_y = TRUE, verbose = TRUE, ...)
{
if (verbose) {
cat(">> preparing promoter regions...\t", format(Sys.time(),
"%Y-%m-%d %X"), "\n")
}
promoter <- getPromoters(TxDb = TxDb, upstream = upstream,
downstream = downstream)
if (verbose) {
cat(">> preparing tag matrix...\t\t", format(Sys.time(),
"%Y-%m-%d %X"), "\n")
}
if (is(peak, "list")) {
tagMatrix <- lapply(peak, getTagMatrix, weightCol = weightCol,
windows = promoter)
}
else {
tagMatrix <- getTagMatrix(peak, weightCol, promoter)
}
if (verbose) {
cat(">> plotting figure...\t\t\t", format(Sys.time(),
"%Y-%m-%d %X"), "\n")
}
if (!(missingArg(conf) || is.na(conf))) {
p <- plotAvgProf.internal(tagMatrix, xlim = c(-upstream,
downstream), xlab = xlab, ylab = ylab, conf = conf,
facet = facet, free_y = free_y)
}
else {
p <- plotAvgProf.internal(tagMatrix, xlim = c(-upstream,
downstream), xlab = xlab, ylab = ylab, facet = facet,
free_y = free_y)
}
return(p)
}

By default, it gives blue n brown colour for 2 samples . I need to change that.
Please help me out. Thank you in advance :)

Answer: Please consider this alternative :

https://github.com/GuangchuangYu/ChIPseeker/issues/89

All the best!
Question: Dear list,

I'm trying to cluster data and then plotting a heatmap using heatmap.2 from ggplot, the script works perfectly with matrices up to 30000 rows, the problems is that I'm using matrices up to 500000 rows (data_sel), and when I try to cluster I get this error:

heatmap.2(as.matrix(data_sel),col=greenred(10), trace="none",cexRow=0.3, cexCol=0.3, ColSideColors=fenot.colour, margins=c(20,1), labCol="", labRow="",distfun=function(x) dist(x,method="manhattan"))
Error in vector("double", length) : vector size specified is too large

Is there any approximation using R to plot heatmpas with his big data?

Thanks in advance

Answer: It simply doesn't work, unless you have around ~930GB+overhead ( almost 1 Terabyte!) of main (free addressable) RAM to hold your distance matrix! Maybe twice or half that much depending on storage mode.
Question: I made a heatmap with the following data (the actual table has 240 rows):

tf.log

0-dpi 1-dpi 7-dpi 14-dpi 22-dpi
TRINITY_DN0_c0_g1_i2 1.27584408 0.5872404 1.780178 1.414362 1.53341851
TRINITY_DN10214_c0_g1_i2 -2.34774378 -2.9194079 -3.211677 -2.885869 -2.04617227
TRINITY_DN10214_c0_g1_i6 -2.14867876 -1.5794957 -1.577811 -2.485095 -1.44172768
TRINITY_DN1038_c0_g1_i4 0.03163921 0.7375222 2.037936 2.462830 0.04559793
TRINITY_DN10462_c0_g1_i2 -2.63973533 -2.8039350 -2.481144 -2.698932 -1.76284020
TRINITY_DN1052_c0_g1_i3 -3.32767605 -2.2006082 -1.887869 -1.211592 -1.43292669

I then made the following heatmap:

heat.gen.k = Heatmap(tf.log, width = unit(10, "cm"),
km = 12,
cluster_columns = F,
show_row_names = F,
row_title_rot = 0,
row_gap = unit(3, "mm"),
name = "Log2FC",
column_title = "Resistant - Susceptible",
column_title_gp = gpar(fontfamily = "sans", fontsize = 28),
column_names_gp = gpar(fontfamily = "sans", fontsize = 20),
col = heat.col,
right_annotation = ha)
print(heat.gen.k)

I would like to pull the clusters out for downstream analysis but I cannot do so with the following code:

clusterlist = row_order(heat.gen.k)
clu_df <- lapply(names(clusterlist), function(i){
out <- data.frame(GeneID = rownames(tf.log[clusterlist[[i]],]),
Cluster = paste0("cluster", i),
stringsAsFactors = FALSE)
return(out)
}) %>%
do.call(rbind, .)

I keep getting the following error:

Error in data.frame(GeneID = rownames(tf.log[clusterlist[[i]], ]), Cluster = paste0("cluster", :
arguments imply differing number of rows: 0, 1

Not sure how to resolve this issue.

Answer: It happens when a cluster consists of only one gene. In that case you are extracting a data.frame with a single row and this is then coerced into a vector which has no rownames. You can easily avoid that by changing:

`rownames(tf.log[clusterlist[[i]],])` into `rownames(tf.log)[clusterlist[[i]]]` so you don't extract the data.frame and then the rownames but extract the rownames first and then subset it to those genes you have per cluster.
Question: Hi,

I want to compare the peaks from two separate chip seq samples.

I'm trying to use MAnorm (R version) on a MacOS unix command line, but I'm having trouble.

Now the program runs through "step1" without any errors, but it doesn't seem to be making a pair of essential files: `peak1.bed` and `peak2.bed`, so all downstream steps fail.

Below, please see the head of each of my test files, and the command I ran.

If you've used this tool before, do you have any idea why I'm having problems?

Thanks very much,

```
$ head t.*
==> t.TreatedPeaks.bed <==
1 1360 3306
1 13366 13931
1 16606 16934
1 20161 21424
1 59752 60275
1 63673 64105
1 67355 67885
1 68243 68654
1 107683 109178
1 116034 117502

==> t.controlPeaks.bed <==
1 20618 21401
1 37881 38291
1 54003 54374
1 57362 59025
1 67346 67844
1 75598 76127
1 76181 76562
1 110638 111047
1 112394 112668
1 117123 117842

==> t.controlReads.bed <==
1 2 52 -
1 11 61 -
1 32 82 -
1 32 82 -
1 32 82 -
1 32 82 -
1 37 87 +
1 37 87 -
1 47 97 -
1 62 112 -

==> t.treatedReads.bed <==
1 1 51 +
1 13 63 +
1 28 78 -
1 31 81 -
1 44 94 +
1 46 96 +
1 64 114 +
1 69 119 +
1 70 120 -
1 77 127 +

$ MAnorm.sh t.TreatedPeaks.bed t.controlPeaks.bed t.treatedReads.bed t.controlReads.bed 266 275
StepI: clean input
StepII: classify common or unique peaks
Error: Unable to open file peak1.bed. Exiting.
Error: Unable to open file peak2.bed. Exiting.
Error: Unable to open file peak1.bed. Exiting.
Error: Unable to open file peak2.bed. Exiting.
StepIII: count peak read
Error: The requested file (read1.bed) could not be opened. Error message: (No such file or directory). Exiting!
```

Here is the usage section of the readme.txt that comes with the MAnorm:

> Usage:

> I. Create a folder and place in the folder MAnorm.sh, MAnorm.r, and all 4 bed files (peak file and read file for the 2 samples under comparison) to be analyzed.

> II. run command: ./MAnorm.sh sample1_peakfile[BED] sample2_peakfile[BED] sample1_readfile[BED] sample2_readfile[BED] sample1_readshift_lentgh[INT] sample2_readshift_length[INT]

> ExampleCommand:

> ./MAnorm.sh sample1_peaks.bed sample2_peaks.bed sample1_read.bed sample2_read.bed 150 150

> The first 4 parameters should be input files in bed format with no header lines
> The first 2 files have ONLY 3 columns: chromosome, start, end.
> The next 2 files should have 4 columns: chromosome, start, end, strand (+/-)

> The last 2 parameters are the number of bp to be shifted for each read. These two parameters are found from MACS peak file *_peaks.xls after "# d =".

Answer: If you look at the shell script, you'll see that it's expecting UCSC chromosome names (i.e., they need to start with "chr"). So just quickly fix that and things should work.
Question: Hi everyone:

I am trying to plot the linkage disequilibrium decay across the genome using genome-wide snps and get a plot that looks something like this:

http://postimg.org/image/5bkdju6jx/

In the methods sections the authors specified that they used the --r2 command in plink. I read the details of this command in the plink website, but I am still not sure how should I use the output file to plot a similar thing.

Any help would be really appreciated.

Thank you in advance.

Javier

Ps: Heres the article's title: Population Genetic Structure of the People of Qatar (2010)

Answer: In the plink.ld file you have colums CHR_A, BP_A, SNP_A, CHR_B, BP_B, SNP_B and R2.

in R do

# install and load ggplot2
install.packages("ggplot2")
library(ggplot2)

# import the data
ld <- read.table("plink.ld", sep="\t",header=T)

# plot the average correlation for each snp distance
ggplot(ld) +
geom_line(aes(x=BP_B - BP_A, y = R2))
Question: Hi,

I am trying to trim the data with trimommatic. Does anyone know how does it work. The manual is not very clear to me. Also I am not clear what adapter file to use, How to set simple mode rather than palindromic mode and what to do with the unpaired files.

I tried to run this command. But I am not sure of what the tool is doing:

java -jar /home/devgen/bin/trimmomatic-0.36.jar PE -phred33 C95VLANXX-2046D-01-01-01_L003_R1.fastq C95VLANXX-2046D-01-01-01_L003_R2.fastq N6_F_P.fastq N6_F_U.fastq N6_R_P.fastq N6_R_U.fastq ILLUMINACLIP:adapters/TruSeq3-PE.fa:2:30:10 MINLEN:25 AVGQUAL:20 -trimlog logs

Thanks for your help.

Answer: Trimmomatic combines a bunch of regularly performed quality control steps in one go.

In your example command, you give it a paired end fastq file with phred33 quality encodings, then clip for the adapters in TruSeq3-PE.fa with a maximum of 2 mismatches with an overlap score of 30 or 10 for single reads. Then, all reads shorter than 25 are discarded (MINLEN:25), and drop all reads where the average quality is below 20 (AVGQUAL:20).

Personally, I prefer a sliding window approach - I often see reads where the quality scores worsen towards the end, it's rare to have a read where the average quality for all bases is below 20, normally the first 2/3 are pretty good. Using a sliding window approach cuts off bad ends.

I'd use (after making sure the reads are actually PHRED33 quality encoded)

PE -phred33 C95VLANXX-2046D-01-01-01_L003_R1.fastq C95VLANXX-2046D-01-01-01_L003_R2.fastq N6_F_P.fastq N6_F_U.fastq N6_R_P.fastq N6_R_U.fastq ILLUMINACLIP:adapters/TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:25

The order of arguments is important - MINLEN25 before the SLIDINGWINDOW approach would mean that it would only check the length of reads before trimming, MINLEN25 after the SLIDINGWINDOW approach means that it will discard reads that after trimming are too short.

If you don't know your adapters, you can run FASTQC on your library - overrepresented k-mers can be adaptersequences. I've also just added all adapters in Trimmomatic's folder of adapters into one file and ran it with that one. There is no distinction between palindromic mode and simple mode - palindromic mode is used for paired end reads, simple for single end reads.

You can use unpaired reads depending on your task - I'd use them for genome assembly since more data is better, but I wouldn't use them for SNP-calling as an alignment of a single read is less precise than an alignment of a read pair.
Question: Hello everyone,

I got a few flanking sequences around a bunch of SNPs. Since they were export from biomart, the variant allele was like a _:

>2|28804|28804|rs544393009|dbSNP|C/T
TGCAT_TCAGC
>2|29358|29358|rs943169482|dbSNP|G/GAAAA
CCCAA_AAGCC
>2|28396|28396|rs781239227|dbSNP|T/C
AATGC_CTTGG

Can some experts help me to replace the _ in each fasta sequence with its header description? That is each fasta input will be turned in 2 or 3 flanking sequence:

>2|28804|28804|rs544393009|dbSNP|C
TGCATCTCAGC
>2|28804|28804|rs544393009|dbSNP|T
TGCATTTCAGC
And that is what I am desperately seeking for!!!

Thanks in advance!

Answer: This was somewhat difficult and posed a challenge. Originally I started with just `sed` and `awk` but it got very complex with them.

Nevertheless, if the format of your input file is *VERY* consistent, the following should work (assuming your input is in *snps.txt*):

paste <(cut -f1 -d_ snps.txt | sed 's/\/[ATGC]* / /g') <(sed 's/[0-9XY]*|[0-9]*|[0-9]*|rs[0-9]*|dbSNP|//g' snps.txt | sed 's/\/[ATGC]* [ATGC]*_[ATGC]*//g') <(cut -f2 -d_ snps.txt) -d ""

2|28804|28804|rs544393009|dbSNP|C TGCATCTCAGC
2|29358|29358|rs943169482|dbSNP|G CCCAAGAAGCC
2|28396|28396|rs781239227|dbSNP|T AATGCTCTTGG

...or, to break it down a bit:

sed 's/[0-9XY]*|[0-9]*|[0-9]*|rs[0-9]*|dbSNP|//g' snps.txt | sed 's/\/[ATGC]* [ATGC]*_[ATGC]*//g' > RefAlleles.list
paste <(cut -f1 -d_ snps.txt | sed 's/\/[ATGC]* / /g') <(cat RefAlleles.list) <(cut -f2 -d_ snps.txt) -d ""