cannin/biostars_qa huggingface dataset:
Question: When running batch file in IGV, how to set its oprion "Show all bases" for track?
Or alternatively, is there an option to show all bases in tracks by default?

Answer: **not tested.** In https://github.com/igvteam/igv/blob/master/src/main/java/org/broad/igv/batch/CommandExecutor.java there is a command: "preference"

https://github.com/igvteam/igv/blob/master/src/main/java/org/broad/igv/batch/CommandExecutor.java#L173

which calls the PreferencesManager

There is a constant **SAM.SHOW_ALL_BASES** https://github.com/igvteam/igv/blob/da21278a6d8980e80a246c7934b549ffb407597c/src/main/java/org/broad/igv/prefs/Constants.java#L123

so I would try:

preference SAM.SHOW_ALL_BASES true

or

preference SAM.SHOW_ALL_BASES 1

?
Question: Hi,

I'm trying to plot a correlation matrix such that the color of each square correspond to the correlation coefficient (deeper color -> higher correlation), while the size of the square corresponds to the pval (larger square -> more significant pval).

Is there any R / python package that can plot such a correlation matrix (as shown)? (The figure is from https://www.nature.com/articles/ng.3406#Sec20)

![enter image description here][1]

[1]: /media/images/78351900-46f9-4e9a-ab7f-0d0aa2ec

Answer: That's `corrplot`: https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
Question: Hello

Simple question ( I hope!)

I'm using this [VariantAnnotation][1] package from bioconductor to work with SNPs.

myparam = ScanVcfParam( # this function is a fast positional filter
samples="3441_pseudo_gatk",
which= GRanges(chr[14], # This is a chromossome name
ranges=IRanges( # these are custom filters
start= 10,
end=10000)))

VCF <- readVcf(TabixFile("my_potato.vcf.gz"), "Potato", param=myparam)

VCF is a "Formal class CollapsedVCF". It's very complex and there are extensive [instructions][2], but I cannot find the answer to the most simple thing of them all:

> How do I parse the raw data of VCF in R?

I know that if I write it out with:

writeVcf(obj = VCF, filename = "tmp.vcf")
I get exactly what I want in a file in my working directory. However, I would like parsing all lines in R as a vector to filter, as **I already wrote a filter that gives index positions**.

All of this because I don't like the filterVcf() function. It doesn't take advantage of the fact that VCF is already filtered by positions and is thus smaller than the original "my_potato.vcf.gz".

**EDIT** : Actually, I got a detail wrong. FilterVcf() can indeed use the positional filter myparam with param=myparam.

Thank you for the attention!
Ricardo

[1]: https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html
[2]: https://bioconductor.org/packages/release/bioc/vignettes/VariantAnnotation/inst/doc/VariantAnnotation.pdf

Answer: > How do I parse the raw data of VCF in R?

You already know the answer to that: it's the `readVcf()` function. You use the `VariantAnnotation` library to fully parse the file, or you use a text file parser if you need the raw input lines for some reason (you very rarely do).

See [here][1] for an example of where I needed to do exactly that to work around a bug in the `VariantAnnotation` library.

> I have a filtering function that outputs desired indexes (based on, for example, the value of info(VCF)$MQ ) but I can't subscript VCF with indexes because it's not a dataframe or a vector.

> I would simply like to filter out certain lines of the vcf. For example, to filter out lines 1,5,7:

Sub-setting a VCF using VariantAnnotation works exactly like you describe. Just subset the VCF using the notation you would use for a vector: `filteredVcf = vcf[-c(1, 5, 7)]` will removes the first, firth, and seventh variant from the VCF.

You can also subset using a logical vector. For example: ` highqualvcf = vcf[VariantAnnotation::fixed(vcf)$QUAL > 20]`

[1]: https://github.com/hartwigmedical/scripts/blob/master/gridss/libgridss.R#L742
Question: I would like to better utilize the functionality available through Bio Perl. I have not been able to find a list of available keywords to use to access different data when using something like Bio::SeqIO. Does anyone have a list like this or can anyone point me to where I would find it? For example: I have written a script to work with fasta files. I needed to check the GI numbers. It was easy to find the keyword accession but this didn't give me the number that I needed. After much online searching I was able to find that `display_id` will give me the number I need. I am sure there are other keywords like `display_id`. What are they?

Answer: With modern Perl object systems, introspection of classes is really easy. You don't have to depend on the documentation (which is usually out of date, and maybe nonexistent), you can use the meta-object methods (the MOP) to introspect the class. Anyway, I'll show you how to inspect classes like those with BioPerl.

The Bio::SeqIO class you mention is really just for IO, or reading and writing, so you probably aren't interested in those methods. You probably want to know what you can do with Bio::Seq objects. Ideally, there would be documentation that lists all of the methods, but BioPerl is a large and mature project, and sometimes you find inconsistencies. In this case, [Bio::Seq][1] is very well documented (just type: "perldoc Bio::Seq"), though I'll demonstrate how to do this anyway. Here is an example to get at what is available when you load the class:

```
use Bio::Seq;
no strict 'refs';

print join ', ', grep { defined &{"Bio::Seq::$_"} } keys %{Bio::Seq::};
```

The hash `%{Bio::Seq::}` is the symbol table of the Bio::Seq class and the `grep` filters out non-subroutine entries (e.g. package variables). This code will tell you the class methods, but I don't think it will give you inherited methods (the code above was modified from [a SO answer][2]).

Here is the result:

```
$ perl biostars108123.pl

validate_seq, primary_id, new, all_SeqFeatures, add_SeqFeature, get_num_of_annotations, DESTROY, feature_count, display_name, remove_Annotations, flush_SeqFeatures, remove_SeqFeatures, desc, description, alphabet, add_Annotation, num_Annotations, primary_seq, top_SeqFeatures, is_circular, namespace, seq, flush_SeqFeature, annotation, accession, id, display_id, object_id, accession_number, subseq, get_SeqFeatures, can_call_new, authority, length, get_Annotations, species, version
```

Not all of those will be of interest to the general public, but most are. If you are unsure about whether your object can do "method_a" you can always test it (with eval or try/catch):

```
eval { $object->can('method_a') };

if ($@) { print "Nope. Can not do 'method_a'"; }
```

[1]: http://search.cpan.org/~cjfields/BioPerl-1.6.923/Bio/Seq.pm
[2]: http://stackoverflow.com/a/1021826/1543853
Question: Hi all,
I have two objectives for my analysis:

1. To find out denovo and rare variants from autism trios analysis `Trios: (Father, Mother, Proband)`
2. To detect frequently occuring mutation across autism patient.

I have used GATK best practice to identify the germline variants across 15 trios (`Total 45 samples [15X3]`) and I have generated `gVCF` file for all the 45 samples.

Can some one guide me how to call `de novo variants` from this trios cohort?

Answer: run GATK VariantAnnotator with '-A PossibleDeNovo --pedigree the.ped '
Question: I have downloaded the complete nr datadase to my local server. I have now a list of nr.00* to nr.38*

nr.00.phd
nr.00.phi
nr.00.phr
nr.00.pin
nr.00.pog
nr.00.ppd
nr.00.ppi
nr.00.psq
nr.00.tar.gz.md5
nr.01.phd
...
nr.38.ppd
nr.38.ppi
nr.38.psq
nr.38.tar.gz.md5
nr.pal
nr.pdb
nr.pos
nr.pot
nr.ptf
nr.pto
taxdb.btd
taxdb.bti

I would like now to run `magicblast` against this local data base, but I keep getting the error massage, I don't have an index file.This is the command I am trying to run:

$magicblast -query R2_val_2.tagged_filter.fastq.gz -db /home/BLASTdb/nr -infmt fastq -outfmt tabular -num_threads 5 -out WK56.R2.unmapped.magicblast.nr.txt

I am getting this error

BLAST Database error: No alias or index file found for nucleotide database [/home/BLASTdb/nr] in search path [/local/projects/Chad/P286::]

Can someone tell me what I am doing wrong?
Do I need to somehow combine these files?

thanks

Answer: yes it does.

my command was as such

magicblast -query tagged_filter.fastq.gz -db BLASTdb/nt/nt -infmt fastq -outfmt tabular -num_threads 5 -out unmapped.magicblast.nt.txt

while the folder with the data base is as such

$ ls -1 BLASTdb/nt
nt.00.nhd
nt.00.nhi
nt.00.nhr
nt.00.nin
nt.00.nnd
nt.00.nni
nt.00.nog
nt.00.nsq
nt.00.tar.gz.md5
nt.01.nhd
nt.01.nhi
...
nt.27.tar.gz.md5
nt.nal
nt.ndb
nt.nos
nt.not
nt.ntf
nt.nto
taxdb.btd
taxdb.bti
Question: I was looking at genome.ucsc.edu for chromosome 1 and realized that it was marked as "NNNNNN" sequences for the end of chromosome.

Why isn't it repeating as "TTAGGG" for these telomeric regions? Isn't N defined as undefined or error regions?

Thank you

Answer: See: https://www.biostars.org/p/66688/#66689

Truly complete sequence of human genome was just pre-published last month. https://www.biorxiv.org/content/10.1101/2021.05.26.445798v1 It will take a while before official genome build gets to that level of resolution.
Question: Hello,
I am new to bioinformatics, and I am trying to do DEG analysis using DESeq2.
I followed the DESeq2 vignette, and I did this

`res <- results(dds, lfcThreshold = log2(2), alpha = 0.05)`

summary(res)

out of 17813 with nonzero total read count.
adjusted p-value < 0.05
LFC > 1.00 (up) : 78, 0.44%
LFC < -1.00 (down) : 106, 0.6%
outliers [1] : 0, 0%
low counts [2] : 2073, 12%
(mean count < 7)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

So, my question is, how do I get my up-regulated and down-regulated genes from this output? I am setting my log2fc=1 for the cutoff threshold.
I also read [this ][1] thread, but I still could not find the answer I am looking for.
Thank you very much.

[1]: https://support.bioconductor.org/p/110866/

Answer: You can convert the `res` object to a dataframe by simply calling `res_df <- as.data.frame(res)`. Then apply filtering like these two functions do:

get_upregulated <- function(df){

key <- intersect(rownames(df)[which(df$log2FoldChange>=1)], rownames(df)[which(df$pvalue<=0.05)])

results <- as.data.frame((df)[which(rownames(df) %in% key),])
return(results)
}

get_downregulated <- function(df){

key <- intersect(rownames(df)[which(df$log2FoldChange<=-1)], rownames(df)[which(df$pvalue<=0.05)])

results <- as.data.frame((df)[which(rownames(df) %in% key),])
return(results)
}

`up <- get_upregulated(res_df)`
`down <- get_downregulated(res_df)`
Question: I am developing an R package and I plan to eventually upload it to BioConductor. However, after I went over the guidelines, I realized that it seems to be more restricted comparing with the requirement from CRAN. I followed the book: "[R Package][1]" by Hadley Wickham, but I believe this book is designed to help people submit CRAN packages, not too specifically for BioConductor.

Therefore I wonder if there is any special issue I should pay additional attention here. Especially formats I guess.

Thank you.

[1]: http://r-pkgs.had.co.nz/

Answer: I have a package at Bioconductor. As a developer you should know:

- Make sure your code is 80 characters wide
- Don't use S3. Bioconductor reviewers will look at your source code.
- Go through much tougher review process than CRAN (but nothing compares to submitting a paper)
- You must have a vignette. This is not the case for CRAN.

Bioconductor packages are usually better than CRAN packages on average, and more user-friendly. As a user, a vignette really helps.

Advantages of Bioconductor:

- Your package will be announced on Twitter by core Bioconductor members. This is a __big-plus__ because they have huge number of followers interested in Bioinformatics.
- Your package will be listed on Bioconductor website. This is useful because that's where Bioinformaticians go.
- You have professional community support
- My paper on my package was submitted and accepted by a respectful journal. The reviewers wrote they trusted the software because it was in Bioconductor.

Disadvantages of Bioconductor:

- Release schedule is 6-months, this is very annoying. Although you can modify the release branch.
Question: I have a set 680k SNP for 141 individuals. Can anyone suggest a software to detect selection?

At the moment I am looking at [gpat++](https://github.com/jewmanchue/vcflib/wiki) and [rehh](http://cran.r-project.org/web/packages/rehh/index.html). Both of them require to specify a region, chromosome and position e.g. 1-1000 for gpat++ and chromosome for rehh (only ~1500 SNPs per chromosome in example). Are there any programs which allow to calculate e.g. iHS for a whole chromosome with large number of SNPs?

Answer: There's a new, very user friendly implementation of most EHH test that you might be interested in.

Here's the link: https://github.com/szpiech/selscan

From the manual:

>
> ```
> selscan -- a program to calculate EHH-based scans for positive selection in
>	genomes
>
> Copyright (C) 2014 Zachary A Szpiech
>
> selscan currently implements EHH, iHS, and XP-EHH, and requires phased data.
> It should be run separately for each chromosome and population (or population
> pair for XP-EHH). selscan is 'dumb' with respect ancestral/derived coding and
> simply expects haplotype data to be coded 0/1. Unstandardized iHS scores are
> thus reported as ln(iHH1/iHH0) based on the coding you have provided.
> ```
Question: Hi

As title says, I am trying to find in which gene a specific region belongs to.

I'm using [solgenomics][1] and for example I'm trying to find in which gene is the region **733456**.

So in **Landmark or Region**: field I put `SL2.40:723456..743456` as content while in the Data source I have selected **ITAG2.4 Release: genomic annotations** and press search.

What server returns me is:

**The landmark named *SL2.40:723456..743456* is not recognized. See the help pages for suggestions.**

So, what am I doing wrong ? Is there any other method to look for that info, or we cannot do such a search?

Thank you

[1]: http://solgenomics.net/gb2/gbrowse/ITAG2.4_genomic/

Answer: Hello Frank,

There is a problem with your search string. You are using "SL2.40:50532..53532" which has 2 issues. ITAG2.4 is using the latest tomato genome build 2.50. The chromosome is also missing. An example of a correct search term is "SL2.50ch08:2626000..2648999". In case you are interested in ITAG2.40, the search string will be "SL2.40ch08:2626000..2648999". Please let me know if you have any more problems.

-Surya
Question: Hello,

I am trying to follow along with the network generation/cluster generation part of the WGCNA tutorial: https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-blockwise.pdf

I am getting the same error regardless of whether I use data for the tutorial or my own RNA-Seq data (both are data frames).

The command that throws the error is (this is the version using tutorial data):

bwnet = blockwiseModules(datExpr, power = 6,
+ TOMType = "unsigned", minModuleSize = 30,
+ reassignThreshold = 0, mergeCutHeight = 0.25,
+ numericLabels = TRUE,
+ maxBlockSize = 2000,
+ saveTOMs = TRUE,
+ saveTOMFileBase = "tstTOM",
+ verbose = 3)

The output I get is:

Calculating module eigengenes block-wise from all genes
Flagging genes and samples with too many missing values...
..step 1
Cluster size 3600 broken into 2133 1467
Cluster size 2133 broken into 1221 912
Done cluster 1221
Done cluster 912
Done cluster 2133
Done cluster 1467
....pre-clustering genes to determine blocks..
Projective K-means:
projectiveKMeans: imputing missing data in 'datExpr'.
To reproduce older results, use 'imputeMissing = FALSE'.
Cluster size 3600 broken into 2133 1467
Cluster size 2133 broken into 1221 912
Done cluster 1221
Done cluster 912
Done cluster 2133
Done cluster 1467
..k-means clustering..
..merging smaller clusters...
Block sizes:
gBlocks
1 2
1973 1627
..Working on block 1 .
TOM calculation: adjacency..
..will not use multithreading.
Fraction of slow calculations: 0.426707
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.
..saving TOM for block 1 into file cnv16pTOM-block.1.RData
....clustering..
....detecting modules..
....calculating module eigengenes..
....checking kME in modules..
..removing 6 genes from module 1 because their KME is too low.
..Working on block 2 .
TOM calculation: adjacency..
..will not use multithreading.
Fraction of slow calculations: 0.358721
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.
..saving TOM for block 2 into file cnv16pTOM-block.2.RData
....clustering..
....detecting modules..
....calculating module eigengenes..
....checking kME in modules..
..removing 1 genes from module 1 because their KME is too low.

Error in (function (x, y = NULL, use = "everything", method = c("pearson", :
unused arguments (weights.x = NULL, weights.y = NULL, cosine = FALSE)

As far as I can tell, the WGCNA library is installed correctly. Can anyone identify what might be causing this "Error in (function..." error? Any ideas as to how to start fixing it?

Thank you!

Answer: It looks like some namespace issue. If you close R Studio, restart your computer, and then reinstall WGCNA as the first thing that you do when you re-open R Studio, then I think that that may resolve it.
Question: Hi
I've done PCA for my gene expression data after DEG analysis, I can see that my case and control samples clustered distinctly but **I want to report a p.value for this result**, it worth to notice that I have the coordinate csv file that represents coordinates of samples from different phenotypes (case vs control).

Thanks in advance

Answer: What Nicolas says is true; however, there are indirect ways to derive a P value based on your PCA results.

#Buid a regression model
You can indirectly derive P values in your situation by building a binary logistic regression model using the PC1 (or PC2, PC3, ..., PCX) values to predict case/control status. From this, you should see a strong P value based on what you say.

summary(glm(CaseControl ~ PC1))

#Correlation
For other type of variables, like continuous variables, you can simply run a correlation test between the continuous variable and the PC1 values, and derive a P value from this too (in R, use `cor.test()`)
Question: sorry,

in mirdeep installation I read

**cd randfold2.0**

**edit Makefile e.g. emacs Makefile**

change line with INCLUDE=-I. to

INCLUDE=-I. -Iyour_path_to_squid-1.9g -Lyour_path_to_squid-1.9g

e.g. INCLUDE=-I. -I/home/Pattern/squid-1.9g/ -L/home/Pattern/squid-1.9g/

does this mean I should go to makefile in ranfold folder and find word INCLUDE=-I and change that to INCLUDE=-I. -

I/home/Pattern/squid-1.9g/ -L/home/Pattern/squid-1.9g or I should type this command in terminal???

Answer: It indeed means to edit the makefile with e.g. emacs. vim, nano, gedit, or any text editor you like
(alternatively, you can also use `sed` but probably with caution)

Make sure to perform exactly the edit as described, unless you'll come across errors soon.
Question: I want to limit the discussion to split-reads.

Say I have a split-read with the primary alignment as `chr1:100-200` and the supplementary alignment as `chr1:300-500`.

Is the SV position I report `chr1:200-300` or `chr1-201-299`

Or something else? I know there is a convention for this, and I wouldn't be asking if I haven't tried googling it.

Thanks!

Answer: The SV position you report depends on the format you report in. I strongly recommend VCF as other formats such as BED and BEDPE are ambiguous for exactly the reason you raise. The encoding of variants into VCF records is covered in the Variant Call Format (VCF) specifications document at https://samtools.github.io/hts-specs/VCFv4.3.pdf. Section 5 is the section you will be interested in.

In your particular example, the variant supported by your split read depends on which, if any, of your split read alignments were aligned to the negative strand and the relative read position of the split reads (i.e. which end is soft/hard clipped). Your read may well support a tandem duplication of base chr1:100-500 inclusive if the CIGAR of your primary alignment is 200S100M and your supplementary is 200M100H.
Question: Hello,

I am trying to get haplotypes for one locus on chromosome19 and for that, I need phased data. However, prior to the phasing, I am supposed to split my PLINK files into separated files by chromosome. I found this script online but it doesn't work (see below). Can you please give me advice on how to do this? I am rather a newbie in the data processing so I am sorry if this is a too basic question.

Thank you!

#!/usr/bin/perl

# This script takes as input the base filename of binary pedfiles (*.bed,
# *.bim, *.fam) and a base output filename and splits up a dataset by
# chromosome. Useful for imputing to 1000 genomes.

chomp(my $pwd = `pwd`); my $help = "\nUsage: $0 \n\n"; die $help if @ARGV!=2;

$infile_base=$ARGV[0]; #base filename of inputs $outfile_base=$ARGV[1]; #base filename of outputs $plink_exec="plink
--nonfounders --allow-no-sex --noweb"; $chr=22; #last chromosome to write out

for (1..$chr) {	print "Processing chromosome $_\n";	`$plink_exec
--bfile $infile_base --chr $_ --make-bed --out ${outfile_base}$_;` }

Answer: No need for scripts, simple loop should work, something like:

for chr in {1..23}; do \
plink --bfile myPlink --chr $chr --make-bed --out myPlink_${chr}; \
done
Question: Hello everybody,

I need to obtain, from a vcf file, a tabular file (to view with excel or libreoffice calc) containing SNP\variant positions, and minor frequency allele.

Is it possible?

Answer: Use the VEP's web interface if you're not tied to the command line:

http://www.ensembl.org/Homo_sapiens/Tools/VEP

or http://grch37.ensembl.org/Homo_sapiens/Tools/VEP if your data is on GRCh37.

See http://www.ensembl.org/info/docs/tools/vep/online/input.html#ident for more details.
Question: Hi Everyone,

I have the data from a ChIP-Seq experiment. I used EDD to call the domains. I am interested in knowing out what genes are present just a few bp outside of the domains. I extended each of the peaks in the bed file I got from edd by 250 bp on both ends using bedtools slop. But what I need now is two rows for each feature of length 500bp centered around the ends. For instance,

Input

chr1 1000 2000
chr1 2500 3500

Output

chr1 750 1250
chr1 1750 2250
chr1 2250 2750
chr1 3250 3750

Using bedtools slop to extend and subtracting the old bed file from new file with extended coordinates will give me the 250bp region exactly outside the two ends of the peak. But I want the ends of the peak to be in the centre. How can I do that ?

Thanks

Answer: A bed file is just a table, so you can do this quickly in R or the language of your choice.

shiftPosition <- function(record, size){
first <- c(record[1], as.numeric(record[2]) - size, as.numeric(record[2]) + size)
second <- c(record[1], as.numeric(record[3]) - size, as.numeric(record[3]) + size)
result <- c(first, second)
return(result)
}

bed <- read.table("your_bed_file.bed", header=FALSE, stringsAsFactors=FALSE)

final <- data.frame(matrix(apply(bed, 1, shiftPosition, 250), ncol=3, byrow=TRUE))

write.table(final, file="your_new_bed.bed", quote=FALSE, row.names=FALSE, col.names=FALSE)

This converts

V1 V2 V3
1 chr1 3680371 3681212
2 chr1 4344252 4350391
3 chr1 4351710 4353125
4 chr1 4491516 4493706
5 chr1 4773006 4777948

to

1 chr1 3680121 3680621
2 chr1 3680962 3681462
3 chr1 4344002 4344502
4 chr1 4350141 4350641
5 chr1 4351460 4351960
6 chr1 4352875 4353375
7 chr1 4491266 4491766
8 chr1 4493456 4493956
9 chr1 4772756 4773256
10 chr1 4777698 4778198