cannin/biostars_qa huggingface dataset:
Question: Is there a tool like bedtools shuffle which I can use to randomly shuffle a bam file? or will I have to convert my bam into a bed and then shuffle it? thanks

Answer: One option is to use *bam2bed* in conjunction with *sample*:

$ bam2bed < reads.bam > reads.bed
$ sample -k 1234 reads.bed > random_sample.bed

One advantage of `sample` over GNU `shuf` is that `shuf` loads everything into memory before shuffling, while `sample` uses reservoir sampling on line offsets so that the memory overhead is much, much lower. Memory usage can be an issue for very large input files.
Question: I am trying to make a co-occurrence network graph for my presence/absence data of genes per genomes but am unsure how to go about with it. I'm hoping to end up with something like the first image below,

Where each gene is linked to another gene , considering if they are both present in the same genomes, where possibly a larger circle being used to describe a higher frequency gene. I originally tried using *widyr* and *tidygraph* packages but I am unsure that my data is not compatible (see second image), as it has the BGCs as rows and the individual genomes as columns. I am examining the presence/absence pattern of the gene pair to determine if they represent a coincident relationship; basically if gene *i* and gene *j* are observed together or apart in the input genomes more often than would be expected by chance.

1) Are there any suggestions on what packages/code I could use that would work with my data set, or how I could adapt my data set to work with these packages?

2) Are there any statistical tests that would be also recommended specifically to assure that there is a coincident or not type relationship?

Network Layout

# Example of data set
# rows = genes
# cols = genomes
set.seed(2222)
df <- matrix(sample(c(TRUE, FALSE), 50, replace = TRUE), 5)
colnames(df) <- letters[1:10]

Thanks in advanced

Answer: To address question 1, I would suggest to use the R [igraph package][1]. There's an [excellent tutorial here][2].
Starting from a binary matrix A that can be considered as the adjacency matrix of the graph, you can do something like:

library(igraph)
G <- graph_from_adjacency_matrix(A)
plot(G)

Here you have a bipartite graph and your matrix is not square so it is not an adjacency matrix but can be considered an incidence matrix. You can expand it to a full adjacency matrix and use the above or you can do:

G <- graph_from_incidence_matrix(A)
Then you just need to style the graph to your liking.

EDIT: Re-reading the question, I see you mean co-occurrence in question 2. There are a number of R packages from different fields that can do co-occurrence analysis from binary matrices such as:
[EcoSimR][3] from ecology (see the [co-occurrence analysis vignette][4]) or [quanteda][5] from text analysis ([tutorial][6]).

[1]: https://igraph.org/r/
[2]: https://kateto.net/network-visualization
[3]: https://CRAN.R-project.org/package=EcoSimR
[4]: https://cran.r-project.org/web/packages/EcoSimR/vignettes/CoOccurrenceVignette.html
[5]: https://CRAN.R-project.org/package=quanteda
[6]: https://tm4ss.github.io/docs/Tutorial_5_Co-occurrence.html
Question: i am using following command to align the reads in hisat2: for my project

[memona@farooq hisat2-2.1.0]$ ./hisat2 –p 64 --max-intronlen 10000 –x /data/memona/hisat2-2.1.0/hisat_index -1 /data/memona/SRR959590_A_1P.fq -2 /data/memona/SRR959590_A_2P.fq –S /data/memona/hisatresult/hisat_align.sam

and getting this error.:

Warning: Output file '64' was specified without -S. This will not work in future HISAT 2 versions. Please use -S instead.
Extra parameter(s) specified: "–x", "/data/memona/hisat2-2.1.0/hisat_index", "–S", "/data/memona/hisatresult/hisat_align.sam"
Note that if files are specified using -1/-2, a file cannot
also be specified. Please run bowtie separately for mates and singles.
Error: Encountered internal HISAT2 exception (#1)
Command: /data/memona/hisat2-2.1.0/hisat2-align-s --wrapper basic-0 --max-intronlen 10000 -1 /data/memona/SRR959590_A_1P.fq -2 /data/memona/SRR959590_A_2P.fq –p 64 –x /data/memona/hisat2-2.1.0/hisat_index –S /data/memona/hisatresult/hisat_align.sam
(ERR): hisat2-align exited with value 1

-p 64 is not a file but number of threads im using. further i want all the output files in hisatresult directory as is spesified in command line..
kindly help me to resolve the issue.

Answer: The `-` in `-p` is an em-dash, not a hyphen, in your command. I imagine you copied and pasted from Word or something like that that "autocorrected" it for you. Just type the command manually and that should resolve the error.
Question: Hi all,

I have few questions about MaSuRCA assembly program and I will greatly appreciate if anyone can clarify these doubts:

Do I have to run trimmomatic or any other filtering/trimming tools before I run the assembly? I am asking this because the program already has a inbuilt error correction step.

How to specify the overlapping libraries ? The sr-config file asks to provide input library as 1)two_letter_prefix 2)mean 3)stdev 4)fastq 1 5) fastq 2, but negative mean are to specify paired-end reads that are outies (RF) and it won't accept 0 as insert size. Also, is the mean, insert-size or the total size (read length * 2 + insert size)?

Do I have to combine similar library types into a single file before putting them in sr-config file or is it OK to specify each library separately?

Thanks very much,

Answer: 1. No need to run trimmomatic. the inbuilt error correction step has adapter trimming, error correction, etc.
2. Read [the quick start guide][1] e.g.

> PE = aa 180 20 /data/fwd_reads.fastq /data/rev_reads.fastq*
The 'mean' and 'stdev' parameters are the library insert average length and standard deviation.

> "JUMP are assumed to be outties <---.--->. If there are any jump libraries that are innies, such as longjump, specify them as JUMP and specify NEGATIVE mean."

So negative mean is for longjump library, not for paired-end reads.

3. No, you can specify each library per line, e.g. add a new line, supposed 250bp PE library with 30bp stdev

PE = bb 250 30 /data/fwd_reads_1.fastq /data/rev_reads_1.fastq

Good luck!

[1]: ftp://ftp.genome.umd.edu/pub/MaSuRCA/MaSuRCA_QuickStartGuide.pdf
Question: I have a EMBL file ready for ENA submission but failing on translation table. It should be 5 and added to file but it seems to expect translation table 1. I can't see where any addition could be to make it pick it up? I added XXX so you don't know the species. Error seems to suggest something on first line, any ideas?

ID XXX; XXX; circular; genomic DNA; XXX; XXX; 15307 BP.
XX
AC XXX;
XX
AC * _Mitochondria
XX
PR Project:PRJEB11111;
XX
DT 01-May-2020 (Rel. 133, Created)
XX
DE XXX
XX
KW .
XX
OS XXX
XX
RN [1]
RP 1-15307
RG XXX
RT ;
RL Submitted (01-MAY-2020) to the INSDC.
XX
FH Key Location/Qualifiers
FH
FT source 1..15307
FT /mol_type="genomic DNA"
FT /organelle="mitochondrion"
FT /organism="XXX"
FT gene complement(12989..13705)
FT /locus_tag="cox2"
FT mRNA complement(12989..13705)
FT /locus_tag="cox2"
FT CDS complement(12989..13705)
FT /locus_tag="cox2"
FT /transl_table=5

Error:

ERROR: organism classified. Submitted /transl_table "5" conflicts with translation table "1" recruited from taxonomy. Please check submitted /transl_table, /organelle and /organism for agreement. Contact us if necessary. [ line: 1 of MT3.embl.gz]

Answer: I think it is inconsistent. In the flat file the sequence is called Mitochondria (cf AC line. It is extracted from the fasta file by EMBLmyGFF3). While in your chromosome list file it is called Mt. You should replace Mt by Mitochondria. Then try again to validate
Question: A Manhattan Plot is created in GWAS studies to visualize where SNP positions and there logarithmic p-values.
Can somebody please with the help of a simple numerical example for 2 or 3 chromosome show how this plot is made?

**Update**:

Actually I am confused how data is processed so that each SNP has different p-values . Take for the sake of argument that I have following scenario

Chromosomes=23
Controls=500
Subjects=800
SNP=20000/SNP that is each chromosome has 20000 SNPs

Now my confusion is that how those dots are made . I mean how is it possible that on horizontal axis that is representing a particular SNP on a specific chromosome how could we have multiple dots going up or down? that is my confusion?

regards

Answer: You seem have some issues reading what is said.

- Each chromosome has multiple SNPs
- For each SNP you perform one statistical test
- Each SNP has one position and one p-value
- Each SNP has only one dot in a Manhattan plot
- It looks like they're lots of dots stacked on each other, but that's just because you are testing many many SNPs.
Question: Hi friends
I plotted this box-wisker for TCGA HTSeq data in R.
I want to have harf of them as red and half as blue (control vs treatment groups). or is there any better way for boxplot?
How can I do that?
I just used simple code:

data <- read.table("mydata.csv", header=TRUE, row.names = 1)
library(DESeq2)
mat <- as.matrix(data)
log <- rlog(mat)
boxplot(log)

my plot:
https://ibb.co/xXsHtC0

Answer: maybe this one? A dodged box plot?

I used "ToothGrowth" public data for this demonstration:

data("ToothGrowth")
library(ggplot2)

ggplot(ToothGrowth, aes(x = as.factor(dose), y = len)) +
geom_boxplot(aes(fill = supp), position = position_dodge(0.9)) +
scale_fill_manual(values = c("#09E359", "#E31009")) +
theme_bw()

![enter image description here][1]

Also, a dogged violin plot would be a very informative way to present expression data. Personally, I would rather use a split violin plot for expression data when I have many cases to be plotted. The codes for creating a split violin plot in python can be found on my [GitHub,][2] which I used to plot top DE genes from RNA-seq.

[1]: https://raw.githubusercontent.com/hamidghaedi/file_sharing/main/box_demo.JPG
[2]: https://github.com/hamidghaedi/RNA-seq-differential-expression/blob/master/violin_plot_DE.ipynb
[3]: https://i.stack.imgur.com/tYQWX.png
Question: I have a network of miRNAs and genes that correlate in my experiment (based on expression levels). The edges link each miRNA to the gene that it correlates with. None of these miRNAs and genes have any known direct interactions according several databases of experimentally validated gene miRNA interactions. It is possible that these correlations exist because of genes which were not measured in my experiment E.g. miR-1 interacts with geneX (not measured in my experiment) which interacts with geneY(measured) leading to a correlation between miR-1 and gene-Y in my network.

I want to create a new network that tries to link (e.g. through gene-gene interactions, perhaps even transcription factors) each miRNA with its correlating gene.

I thought cytoscape may be able to help but I have never used it before and am struggling to find a plug-in which does what I want.

I have a list of experimentally validated gene miRNA interactions and I presume it is possible to get a list of gene-gene interactions e.g. from an interactome database. Does anyone have any ideas on how can I integrate this information with my network so that I can link each miRNA with its correlating gene?

Many thanks

Answer: Greetings, and welcome to Cytoscape. I think there are a couple of plugins that might be useful. What I would suggest is to use your list of genes and then probe one of the databases to explore how those genes might be related. For example, you could use the stringApp, add your list of genes, but then also add a number of high confidence edges that might bring in interacting genes that links them. If you want to view your miRNA interactions in the same network, Cytoscape has an excellent "merge networks" tool that will allow you to merge your original network with the STRING network. You might also want to look at the new IntActApp as an alternative using the IntAct database in lieu of the STRING database.

-- scooter
Question: I am preparing admixture analyses involving the Neanderthal and Denisovan genomes, and I have downloaded extended (=generally non-vcftools friendly) VCF files (http://cdna.eva.mpg.de/neandertal/altai/AltaiNeandertal/VCF/ and http://cdna.eva.mpg.de/denisova/VCF/hg19_1000g/). The files were originally made with GATK, but the authors greatly modified to files; thus, they aren't standard VCFs.

I've got them down to just the sites that I have modern data for, so they only a fraction as massive as the original files.

A few hundred sites sites have been marked LowQual, and have been ejected using `grep -v`; however, LowQual is going into more than just Genotyping Quality (i.e., LowQual sites have GQ's as high as 59, but for the Neanderthal 8,767 of 511,858 non-LowQual sites have GQ<60).

I was wondering what would be good GQ cut-off to use for the non-LowQual line?

Also, any suggestions on how to filter them? (Remember I cannot use VCFtools or GATK or anything similar due to the non-standard formatting)

Here is an example line, GQ is in the subsequent code block

1 5031561 rs7518523 A G 909.02 . AC=2;AF=1.00;AN=2;DP=24;Dels=0.00;FS=0.000;HRun=0;HaplotypeScore=0.9665;MQ=39.43;MQ0=0;QD=37.88;1000gALT=G;AF1000g=0.34;AFR_AF=0.70;AMR_AF=0.28;ASN_AF=0.17;EUR_AF=0.27;UR;TS=HPGOMC;TSseq=A,G,A,G,G,A;CAnc=A;GAnc=A;OAnc=G;bSC=987;mSC=0.000;pSC=0.007;GRP=-2.27;Map20=1 GT:DP:GQ:PL:A:C:G:T:IR 1/1:24:72.23:942,72,0:0,0:0,0:10,14:0,0:0

```
GT 1/1
DP 24
GQ 72.23
PL 942,72,0
A 0,0
C 0,0
G 10,14
T 0,0
IR 0
```

Answer: Since I didn't get any suggestions, I sought out how this dataset has been used recently in the literature.

I found sources by Qin and Stoneking (dx.doi.org/10.1093/molbev/msv141) and Lazaridis et al. (dx.doi.org/10.1038/nature13673) (former cites the latter), which suggest filtering the LowQual sites as well GQ < 30 and Qual < 50.

Just thought to pass this along for anyone else using these datasets.
Question: how can i remove duplicated variants from vcf file? i googled and searched in biostars history but i did not fond any way to do it.

Answer: use [vcfuniq][1] or [bcftools][2] norm (with `-d` option) to remove duplicates

[1]: https://github.com/vcflib/vcflib#vcfuniq
[2]: https://samtools.github.io/bcftools/bcftools.html
Question: Hi.

I have a txt file containing mutiple fasta sequences, and I'd like to replace Nth character in each sequence.

If I want to replace 10th character of each sequence in the example below, which linux command can I use?

> sample1
TATCCGATGCGACGTGCAGCG
> sample2
CTAGCGTAGTGTCGACTGCAT
> sample3
GACTGACGTGACGTAGTCGAC

Thank you!

Answer: Assuming one-line sequences in your FASTA-formatted input, replace `X` with your character of interest:

$ awk '{ \
if ($0 ~ /^>/) { \
print $0; \
} \
else { \
printf("%s%c%s\n", substr($0, 1, 9), "X", substr($0, 11, length($0) - 10)); \
} \
}' in.fa > out.fa

If you have multi-line sequences in your FASTA input, you need to render them into single-line sequences before using this one-liner. See the following Biostars question (and answer) for a solution: https://www.biostars.org/p/9262/
Question: I have a bed file and wanted to calculated average PhastCon score over the regions, is there a tool which could do so? I tried using cmotifs but it submits jobs and at times does not produce timely results.

Answer: You could use *bedmap --mean* here, to map PhastCon signal onto your sorted regions.

First, grab PhastCon data from UCSC and convert WIG to a compressed form of BED called Starch, using a tool called *wig2starch*:

```
$ rootPath="http://hgdownload.cse.ucsc.edu/goldenpath/hg19/phastCons100way/hg19.100way.phastCons"
$ for i in `seq 1 22` X Y; \

do \
echo "getting PhastCon data for chromosome chr${i}..."; \
wigFn="chr${i}.phastCons100way.wigFix"; \
url="${rootPath}/${wigFn}.gz"; \
wget -qO- ${url} | gunzip -c - | wig2starch - > ${wigFn}.starch; \
done
```

We could have used `wig2bed` instead of `wig2starch`, but note that the final set of PhastCons files could take between 10-20 GB compressed, and uncompressed BED files could require about 10x more disk space. It takes a little longer to make Starch files, but it will take less disk space. I'd recommend some kind of compression, if you plan to do repeated map calculations.

Second, map PhastCon signal files to your sorted regions file by chromosome:

```
$ for i in `seq 1 22` X Y; \
do \
echo "mapping signal for chromosome chr${i}..."; \
phastConFn="chr${i}.phastCons100way.wigFix.starch"; \
bedmap --echo --mean --chrom chr${i} regions.bed ${phastConFn} > regions_with_avg_phastcons.${i}.bed; \
done
```

Third (optionally), union the per-chromosome results to a single file, using *bedops --everything*:

$ bedops --everything regions_with_avg_phastcons.*.bed > regions_with_avg_phastcons.bed

Doing the *bedmap* step with `--chrom` allows parallelization with grid job schedulers or GNU Parallel, splitting the tasks by chromosome. This could reduce calculation time down to that required for the largest chromosome of regions.

One thing I would note is that the average value calculated is over scores for mapped elements. If you have gaps of scores along your mapped region, those gaps would not contribute to the calculation of the mean signal.

I do not think you can interpolate PhastCon conservation signal between a gap, but it has been a while since I have worked with this particular dataset - its authors could probably confirm or negate this.

It might generally be useful to use *bedmap* with the *--skip-unmapped* option to filter out results for unmapped regions.
Question: Is anyone aware of genome-wide CRISPR/Cas9 screens where gene expression is quantified? Meaning, after targeting a gene, what is the average gene expression change for all genes (microarray, RNAseq, etc.)?

Thanks in advance

Answer: There are some papers, done in specific cell types and disease models:

https://www.nature.com/articles/nature23477?_ga=2.135383345.140789132.1527763165-716992680.1527763162

https://www.cell.com/cell-reports/pdf/S2211-1247(18)30387-5.pdf

https://www.ncbi.nlm.nih.gov/pubmed/29038160
Question: Hello,

Is there a script somewhere that I can use to convert the `FRG `files back to paired end `fastq `files?

If not, what info can I use to find the pairs? First 2 records for my `FRG` file are as follows:

```
{FRG
act:A
acc:aa202
rnd:1
sta:G
lib:aa
pla:0
loc:0
src:
.
seq:
GCACATGTGCATAGATTTCACGGACCTCAATCGCGCCTGTCCGAAAGATGACTTCCCTTTGCCCCGGATAGATCAACTGGTCGATTCTACGGCCGGCTGCGAAGCGATGAGTTTCTTGGATGCTTATTCTGGCTACCACCAGATTAGCAT
.
qlt:
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
.
hps:
.
clr:0,150
}
{FRG
act:A
acc:aa203
rnd:1
sta:G
lib:aa
pla:0
loc:0
src:
.
seq:
TTGCCCTTGAAAGGTCTTTGTCATGAGTCTCTGGAAGGTTGCTCCGGCGTTTTTCAAGCCGAATGGCATTCGGACGTAGTAGAAAGTACCCACAGGAGTGATGAAACTAGTCTTCTCTTCATCTTCCGGGATCATGCTAATCTGGTGGTA
.
qlt:
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
.
hps:
.
clr:0,150
}
```

Thanks for any help!

EDIT: I can easily parse the FRG file to fastq, but still don't know how to find pairs. If you any of you have experience working with frg/celera assembler, please provide any insights. Thanks.

Answer: There is a conversion script in celera:

celera/src/AS_GPK/frg-to-fastq.pl

Judging from the output file names, it is PE aware.

my $outAname = "$outPrefix.1.fastq";
my $outBname = "$outPrefix.2.fastq";
my $outIname = "$outPrefix.i.fastq";
my $outSname = "$outPrefix.s.fastq";

There is also one in the AMOS package (not sure if this does PE though)

amos/bin/frg2fastq
Question: > Several very commonly used annotation databases for human genomes are additionally provided below. In general, users can use -downdb -webfrom annovar in ANNOVAR directly to download these databases. To view of full list of databases (and their size and last changed date) **prepared by ANNOVAR developers, use avdblist keyword in -downdb operation**

Ok so I tried this

perl annotate_variation.pl -downdb avdblist

And this

perl annotate_variation.pl -downdb -webfrom annovar avdblist

And this

perl annotate_variation.pl -downdb -webfrom avdblist

And other variations and I get no list

WTF is it so hard to give an example of the exact command to use?

Answer: Thanks to natasha.sernova I scrolled down a little bit and found the answer (and someone just as frustrated as me)

This works

$ perl annotate_variation.pl -webfrom annovar -downdb avdblist -buildver hg19 .

![ANNOVAR needs better documentation][1]

[1]: http://i.imgur.com/Y5OWXoA.png
Question: I am planning to use [bam2fastx][1] tool for converting the BAM files from PacBio sequencing to fastq. (I prefer a command line tool over SMRT link). However, I cannot find any manual for this tool.

The help says that there is an option called `--split-barcodes` for demultiplexing. However, there is no information about the format to provide the barcodes in. My sequencing data contains ~120 asymmetrically PCR-barcoded samples. How do I split the data into multiple files using this tool?

[1]: https://github.com/pacificbiosciences/bam2fastx/

Answer: Documentation for the commands can be found in the official [SMRTanalysis guide][1]. I don't see any way to specify index sequences so it must be recognizing official PacBio indexes. There is [one more page here][2].

[1]: https://www.pacb.com/wp-content/uploads/SMRT_Tools_Reference_Guide_v600.pdf
[2]: https://devhub.io/repos/PacificBiosciences-bam2fastx
Question: Hi everyone,

I have a number of Human Heart RNASeq samples. I have generated a mpileup file on each of them using a bed file containing a list of positions for which I want the mpileup to be made. I have included the -B option to turn off the BAQ computation in mpileup.

Single sample mpileup:

samtools mpileup -uf hg19.fa -B -l snps_chr17.bed sample1.bam > sample1.mpileup

Multiple sample mpileup (all samples):

samtools mpileup -uf hg19.fa -B -l snps_chr17.bed -b bam_files_list.bam > all_samples.mpileup

Now, I want to output the frequency of nucleotides A, T, G and C at each of those positions. There is a perl script pileup2baseindel which generates exactly the output that I want, only that it is incorrect (inconsistent with what is seen in IGV). This program takes a single/multi-sample mpileup file and produces the following output. But like I said, the output is incorrect.

I have searched a lot but can't seem to find any program that does the same. Does anyone know any tool/script that would give me such an output with/without taking the mpileup file as an input. Suggestions on alternatives/options are welcomed.

[UPDATE] There was a bug in the perl script and I have updated the code, it works fine now. I have contacted the author and as soon as the author responds, I will send it to him so that he can update it. I would still like to get suggestions and try out other programs.

[UPDATE] Thanks everyone who helped me with this problem. Each of your suggestions worked but I accepted Chris Miller's answer recommending the use of bam-readcount because of its easy of use and the elaborate output that it generates. It directly takes, as input, a bam file as well as a list of positions in bed format for which you want the nucleotide frequency. So you don't need to create an intermediate mpileup file - as in pileup2baseindel or run the program against positions that you are not interested in - as in pysamstats.

Thanks!

Answer: Another option is just to run [bam-readcount][1], giving it the bams and specified positions. It's more straightforward than creating an mpileup and parsing that ugly string.

[1]: https://github.com/genome/bam-readcount