cannin/biostars_qa huggingface dataset:
Question: Hello BioStars community,

New to perl and programming in general, so I thought I might try out my luck asking a question here.

I am trying to match a fairly conserved protein sequence to a proteome using a regex. I am able to output the matching lines, as well as their positions, but I cannot find a way to output the accession numbers along with lines that match my conserved protein.

Here's part of my code:

my $proteins;
open( file, "Athaliana_167_protein.fa" ) or die "can't open file!";
while (){
if (/W[S]TRRKIAI/) {print}
}

Would using lookahead/lookbehinds possibly work to print out the match line and accession number?

Thanks!

Answer: This statement is important for parsing fasta efficiently in perl (unless you want to use BioPerl):

local($/) = "\n>";

It allows you to read a complete fasta record instead of each line.

{
local($/) = "\n>"; # read each fasta record, always use local!
while (my $fastarec = ) {
chomp $fastarec;
my ($defline, @seq) = split "\n", $fastarec; #seq id is the first line
$defline =~ s/^\>//; # remove left over >, just in case
my $seq = join "", @seq; # put together the sequence again
$seq =~ s/\s//g; # remove potential left-over spaces, empty lines etc.
if $seq =~ /W[H|A]TEVER/ {
print (join "\n", ">$defline", @seq), "\n"; # output sequence in original formatting
}
}
}
Question: I was wondering if anyone knew of new software that can help find novel genes and isoforms from my RNA-seq data? I know there Cufflinks and what not but that is a bit old. I was wondering if anyone knew of anything newer? A lot of the genes in the genome that I am working with are predicted. I was interested to see if the RNA data matched up with the gene prediction?

Answer: [StringTie][1] is the sucessor of Cufflinks. There are plenty of relatively recent programs, for example, [KNIFE][2], [BIISQ][3], or [SparseIso][4].

P.S.: see a comprehensive list at [Omicstools][5].

[1]: https://ccb.jhu.edu/software/stringtie/
[2]: https://github.com/lindaszabo/KNIFE
[3]: https://github.com/bee-hive/BIISQ
[4]: https://github.com/henryxushi/SparseIso
[5]: https://omictools.com/novel-transcript-quantification-category
Question: I have a VCF file where in the INFO header I see `Number=1` and `Number=0`. What does that mean? That one is printed and not the other?

##contig=
##ALT=
##INFO=
##INFO=
...

Why in the INFO column, I can't see INDEL, IDV, IMF and that it starts at DP?

#CHROM POS ID REF ALT QUAL FILTER INFO ...
scaffold440 420 . T G 999 . DP=175;VDB=4.12831e-08;SGB=14.3961;RPB=0.362442;MQB=0.00113257;MQSB=6.33575e-07;BQB=0.21253;MQ0F=0;AF1=0.07
scaffold440 451 . C A 999 . DP=174;VDB=2.7432e-08;SGB=4.37813;RPB=0.00386186;MQB=0.259306;MQSB=1.3575e-06;BQB=0.237116;MQ0F=0;AF1=0.102

Answer: The VCF format is kind of all over the place. The fields are generally not required. Rather, those lines tell you what the field means, if it happens to be present.

"Number=1,Type=Integer" indicates that if that tag is present it will be followed by a single integer value (e.g. "DP=175"). "Number=0,Type=Flag" means it won't be followed by anything (e.g. "INDEL") because the presence of the tag itself indicates some kind of information (like in this case that the variant call is an indel, in case you are unable to determine that yourself by noticing that the length of the ref and alt are different).

Edit: I would be remiss to not mention the [official VCF spec][1] in this post.

[1]: https://samtools.github.io/hts-specs/VCFv4.3.pdf
Question: Dear all,

I have indexed the C. elegans reference genome with:

bwa index output/genome/ref/seq/celegans.fa

and then aligned my de novo assembly to the reference with:

bwa mem -t 8 -x intractg output/genome/ref/seq/celegans.fa input/assembly/celegans/hgap/bristol/assembly.fa > output/alignment/pacbio/bwa/ref/bristolAssembly.sam

I then

samtools view -bS output/alignment/pacbio/bwa/ref/bristolAssembly.sam > output/alignment/pacbio/bwa/ref/bristolAssembly.bam

samtools sort output/alignment/pacbio/bwa/ref/bristolAssembly.bam -o output/alignment/pacbio/bwa/ref/bristolAssemblySorted.bam

samtools index output/alignment/pacbio/bwa/ref/bristolAssemblySorted.bam

to visualise the alignment in IGV.

This is what I get: [http://cristian-riccio.ch/wp-content/uploads/2017/07/igv_snapshot.png][1]

You can see that my assembly is the same as the reference, just shifted by 12 bases or so. Has anyone any suggestions about how to solve this problem?

Thanks.

Best,

C.

[1]: http://cristian-riccio.ch/wp-content/uploads/2017/07/igv_snapshot.png

Answer: Solved finally!

[http://cristian-riccio.ch/wp-content/uploads/2017/07/igv_snapshotSolution.png][1]

I used this command instead, penalising gap extensions less.

bwa mem -t 8 -E0.5 -x intractg output/genome/ref/seq/celegans.fa input/assembly/celegans/hgap/bristol/assembly.fa > output/alignment/pacbio/bwa/ref/bristolAssembly.sam

The assembly comes from pacbio, and this technology is 12% indel error in the raw reads, that might be why the contigs also need a bit more lee-ways on indels for alignment.

I figured it out with a colleague because we looked at the start of the chromosomes and it started off aligned and got misaligned when there was a larger indel in my assembly.

Thanks to all for the suggestions.

Best,

C.

[1]: http://cristian-riccio.ch/wp-content/uploads/2017/07/igv_snapshotSolution.png
Question: I found some tips for speeding up `read.table()` [here][1], I wondered if anyone can suggest something for [read.fasta()][2], which is part of the seqinr library. Not that they're really comparable.

I have a ~6GB file.

[1]: http://www.biostat.jhsph.edu/~rpeng/docs/R-large-tables.html
[2]: http://www.inside-r.org/packages/cran/seqinr/docs/read.fasta

Answer: Since that function is written in R I don't think optimizing your code will result in much of a speed up compared to using a library written in another language. If you want to work in R, it would probably be easier to use the `system()` command to transform your sequences with another program or just use command line as Pierre suggested. If you have a specific task I'm sure people can provide a solution.
Question: This dataset is for a specific disease-gene-test results. The dataset goes like this.

| Test | Gene | Relevance | Values |

Test and Gene are the two parameters on x and y axis. Values are the combined results for the parameters pair.
The problem for me here is that I have one more parameter called Relevance which represents the relevance of the test-gene pair and it is boolean (only two values YES/NO).

The dataset (Relevance) should be differentiated in the map with different colour (like red and green).

Gradient of that colour represents numerical values between that interaction.

The end result I was aiming at was, Test on x-axis and Gene on y-axis and the map for these interaction will be only with two colours (representing Relevance values) and the gradient of that colour representing Values.
Is this possible to achieve this kind of Heatmap, if Yes how can I achieve.
If not is there any other option to display such kind of data (similar to heatmap)

Help appreciated !!

Thanks,

RDS

something like this -

Answer: Here is some R code that may help. I must admit I do not understand exactly how your Relevance and Value variables are related (or expected to interact). I have made some guesses. If I have guessed wrong, perhaps you can post some sample data to help clear up my confusion?

I have used scale_fill_gradient2 in these examples. You can specify three different colors: high, mid (default is white), low. You can specify the midpoint (value that maps to mid color, default=0), and upper and lower value limits. This may provide enough flexibility to show your data the way you want it.

In Example 1, negative values range from red to white, and positive values range from white to blue. In Example 2, where relevance is FALSE, no color is plotted and where relevance is TRUE, the color gradient spans the range red-white-blue.

library(ggplot2)

# Create test data.
dat1 = data.frame(x=factor(rep(c("A", "B", "C"), 3)),
y=factor(rep(c(37, 8.7, -17.7), c(3, 3, 3))),
z=c(34, 18, 31, 9, -2, 4, -21, -33, -13))

p1 = ggplot(dat1, aes(x=x, y=y, fill=z)) +
theme_bw() +
geom_tile() +
geom_text(aes(label=paste(z))) +
scale_fill_gradient2(midpoint=0, low="#B2182B", high="#2166AC") +
opts(title="Example 1")

ggsave(plot=p1, filename="plot_1.png", height=4.5, width=5)

# Create a slightly different test dataset.
dat2 = data.frame(
gene=factor(rep(c("Gene_A", "Gene_B", "Gene_C"), 3)),
test=factor(rep(c("Test_1", "Test_2", "Test_3"), c(3, 3, 3))),
relevance=c(TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, TRUE),
value=c(-16, 3, NA, NA, -13, 25, -4, NA, -26))

p2 = ggplot(dat2, aes(x=gene, y=test, fill=value)) +
theme_bw() +
geom_tile() +
geom_text(aes(label=paste(value))) +
scale_fill_gradient2(midpoint=0, low="#B2182B", high="#2166AC") +
opts(title="Example 2")

ggsave(plot=p2, filename="plot_2.png", height=4.5, width=5)
Question: There are two VCF files that I like to merge them, using GATK or VCFtools. The problem is, they have different chromosomal notation, one has Chr, the other does not. This question could be similar to [this one](/p/13462/)

Is there any quick awk/sed commands that you suggest ?! Also I appreciate if you make comment, which of these two (GATK/VCFtools) is more reliable for this task.

Answer: awk '{gsub(/^chr/,""); print}' your.vcf > no_chr.vcf

Should get rid of the chr

awk '{if($0 !~ /^#/) print "chr"$0; else print $0}' no_chr.vcf > with_chr.vcf

Will add the chr to the VCF without chr.

You can use either command depending on how the chromosomes are named in your reference.

Regarding preference of tools, if you plan to do downstream processing with GATK, I'd suggest sticking with GATK CombineVariants for consistency.
Question: Duplicate reads have first been removed using picard:

java -jar -Xmx3g picard/dist/picard.jar MarkDuplicates INPUT=input.bam OUTPUT=output.bam METRICS_FILE=output.dup_metrics CREATE_INDEX=TRUE VALIDATION_STRINGENCY=SILENT

When I run samtools flagstat on the output bamfile I get the following:

2182812 + 0 in total (QC-passed reads + QC-failed reads)
226710 + 0 duplicates
2176925 + 0 mapped (99.73%:-nan%)
2182812 + 0 paired in sequencing
1091406 + 0 read1
1091406 + 0 read2
2156992 + 0 properly paired (98.82%:-nan%)
2171322 + 0 with itself and mate mapped
5603 + 0 singletons (0.26%:-nan%)
9776 + 0 with mate mapped to a different chr
7030 + 0 with mate mapped to a different chr (mapQ>=5)

So presumably the file still contains 226710 duplicate reads.
If I filter out duplicates according to this table:

Flag Chr Description
0x0001 p the read is paired in sequencing
0x0002 P the read is mapped in a proper pair
0x0004 u the query sequence itself is unmapped
0x0008 U the mate is unmapped
0x0010 r strand of the query (1 for reverse)
0x0020 R strand of the mate
0x0040 1 the read is the first read in a pair
0x0080 2 the read is the second read in a pair
0x0100 s the alignment is not primary
0x0200 f the read fails platform/vendor quality checks
0x0400 d the read is either a PCR or an optical duplicate

using the command:

samtools view -F 400 output.bam | wc -l

I get 506072,
not: total reads (2182812) - duplicates (226710) = 1956102

My question is why does samtools flagstat indicate that there are still duplicates present after running picard tools, and why are these figures inconsistent when I attempt to filter out duplicates using samtools? I've been asked to remove duplicates for a project. At the moment I am very confused as to which method I should use given the inconsistencies in results.

Answer: >samtools view -F 400 output.bam

To filter out duplicates use `-F 1024` or `-F 0x0400`, not 400. See also https://broadinstitute.github.io/picard/explain-flags.html
Question: Hi,
I have a list of genes from a Pseudomonas strain with official gene identifiers and I have been trying to find a resource in order to sort this genes according to their metabolic pathways. Also if possible I'd like to also represent these results in some sort of graph like a circular one or a barplot but this is secondary. This strain is not supported by databases like KEGG but it is included in the NCBI. I have tried R packages like clusterProfiler and web servers like DAVID 6.8 and I have not been able to make it work out so far.

Any tips/advice would be greatly appreciated.
Cheers
PS: I have recently started learning R and UNIX so I'm not very proficient.

Answer: There are indeed a lot of tools and resources to choose from. Three important considerations are (1) methods, (2) update frequency, and (3) coverage.

1. METHODS. The methods for pathways classification you are asking about are generally refered to as functional enrichment analyses and include over-representation analysis (ORA), gene set enrichment analysis (GSEA), and topological analysis. Each resource implements one or more of these methods usually with a bit of customized pre or post filtering. ClusterProfiler, for example, supports both ORA and GSEA. There are review papers that compare these methods, but as there aren't really any gold standard benchmarks, many folks just try more than one method. In any case, it's important to know which method you are using with a give tool and its caveats.

2. UPDATES. In terms of sources of pathway annotations: KEGG has not been significantly updated since around 2011 (https://www.kegg.jp/kegg/docs/upd_map.html), BioCyc is updated around 3 times a year (https://biocyc.org/release-notes.shtml), and WikiPathways is significantly updated every month (http://releases.wikipathways.org)! So, this is something to definitely take into account when choosing a resource to plug into a given method. In terms of tools: DAVID hasn't been updated since 2016 (https://david.ncifcrf.gov/gene2gene.jsp), which means their GO and pathway annotations will all be outdated. ClusterProfiler is under active development (https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) and let's you supply your own GMT file from WikiPathways or wherever to perform ORA or GSEA.

3. COVERAGE. If you're studying human or mouse models, then pretty much any resource will do. But you are working with a Pseudomonas strain. This will dramatically limit your options. In general, I would suggest BioCyc for baterial species coverage. It looks like they currently support two different strains for Pseudomonas. You might consider translating your genes to E.coli identifiers to perhaps expand your options to other resources.
Question: I would like to use the Gviz package but my I have a reference genome in FASTA format, not as a valid UCSC genome (it was not published yet).

It is possible to use this FASTA file instead?

Let´s say that I have one genome in FASTA format on disk:

library(Biostrings)
ncrna <- readDNAStringSet(file = "GTgenome.fa")

head(ncrna)
A DNAStringSet instance of length 6
width seq names
[1] 20202851 CCCAGTTTTCCCCACTCTGTGA...AAAGATCTTACAACCGATTTT chr10 Major...
[2] 20315886 AGCCGACGAGACTCACAGAACC...TCACAAACCCCCTCGGGAGGG chr11 Major...
[3] 20466350 GATTAGACCTCCGAAAGGGGTA...ATTATTAATTATTAAATATTA chr12 Major...
[4] 16480340 GTCTCCACTTGCCCCACAACGG...AGATGACGATGATGAAGATGA chr13 Major...
[5] 16193477 CTCTGTGACATCACAGCCATGG...GGGTTACACACGTTGTTTTTT chr14 Major...

Using the object created to replace a UCSC genome:

ideoTrack <- IdeogramTrack(genome=ncrna, chromosome="chr10", fontsize=14)

Error in .Call2("new_XStringSet_from_CHARACTER", ans_class, ans_elementType, :
key 51 (char '3') not in lookup table
In addition: Warning message:
In if (!token %in% base::ls(env)) { :
the condition has length > 1 and only the first element will be used

I am probably misunderstanding a very basic feature here, but I would be grateful for some help!

Answer: I am not sure if I fully understood the problem, but what you could try is to replace an UCSC Ideogram with your own data. This solution is certainly not the best way, but I'm not that familiar with the Gviz package.

Build a track from, for example the human genome chr22. Subsequently, replace the necessary slots (range, chromosome bands and names) with your own data.

library(Gviz)
ideoTrack <- IdeogramTrack(genome = "hg38",chromosome = "chr22")
plotTracks(ideoTrack,from = 1,to = 1000000)

# set your own
my_ideoTrack <- ideoTrack

# RANGE
START=seq(1,100,10)
END=seq(10,100,10)
RANGES <- GRanges(seqnames = c(letters[1:10]),IRanges(start =START ,end =END))

# arbitrary Giemsa staining
GS <- c("gneg", "gpos100", "gpos25","acen", "gpos50", "gpos75", "gvar", "stalk" ,NA,NA)
bt <- data.frame(chrom="own1",chromStart= seq(1,100,10),chromEnd=seq(10,100,10),name=letters[1:10],gieStain=GS)

# chromosome name
CHR <- "own1"

# Replace with your own data
my_ideoTrack@range <- RANGES
my_ideoTrack@chromosome <- CHR
my_ideoTrack@name <- CHR
my_ideoTrack@bandTable <- bt

# set UCSC chromosome Names False
options(ucscChromosomeNames=FALSE)
gx <- GenomeAxisTrack()
# Plot
plotTracks(list(my_ideoTrack,gx),from = 1,to = 15)
Question: I have analyzed RNA-seq data with DESeq2 and am trying to plot a 3D PCA using rgl-plot3d.

I was trying to output PC1, PC2, and PC3 and then plot them. However, I realized that I get different results for PC1 (and PC2) when I try plotPCA (used with DESeq2) and prcomp.

What is the bug on my code?

dds <- DESeqDataSetFromHTSeqCount( sampleTable = sampleTable,
directory = directory,
design= ~group)
rld <- rlog(dds, blind=TRUE)

**From DESeq2:**

data <- plotPCA(rld, intgroup=c("treatment", "sex"), returnData=TRUE )
data$PC1
> [1] -1.9169863 -2.0420236 -1.9979900 -1.8891056 0.9242008 1.0638140

>[7] 0.6911183 1.0551864 0.9598643 -1.5947907 -1.5666862 -1.6694684

>[13] -1.2523658 -1.0785239 1.3005578 2.2913536 2.5381586 2.4287372

>[19] 1.7549495

**Using prcomp**

mat <- assay(rld)
pca<-prcomp(t(mat))
pca <- as.data.frame(pca$x)
pca$PC1

>[1] -1.29133735 -2.96001734 -3.08855648 -3.51855030 -0.68814370 -0.01753268

>[7] -2.31119461 -0.10533404 -1.45742308 -1.30239486 -1.36344946 -1.93761580

>[13] 6.04484324 4.83113873 0.75050886 -0.14905189 2.70759465 3.43851631

>[19] 2.41799979

Answer: Before calling `prcomp()` internally (HERE), the DESeq2 function pre-filters your data based on variance. There is a parameter that can be modified to indirectly disable this behaviour..

Kevin
Question: Hey, I am actually very new to this and I have few genomes to analyse. I am looking through many methods, could anyone please suggest me any good method (reliable but relatively easy to use) for identifying repetitive DNA/ transposons?

Answer: If you have large sequences to analyse, you can download the RepeatMasker app and it will be able to identify the different classes of transposons and repeat elements. RepeatMasker requires a proprietary sequence database called repbase, which can be restrictive. If you want to avoid repbase/repeatmasker there are now some alternatives like DFAM. Links below.

https://www.repeatmasker.org/RepeatMasker/

https://github.com/HullUni-bioinformatics/TE-search-tools

https://dfam.org/home

https://bioinformatics.stackexchange.com/questions/343/are-there-any-repbase-alternatives-for-genome-wide-repeat-element-annotations
Question: Hi everybody,
I have some CNV data that need to be turned into a chromosomal plot. The format is like this:

chr start end ploidy loss/gain

chr1 86000000 117150000 1 loss

chr2 70250000 70500000 3 gain

chr2 203050000 204650000 3 gain

(The last column can probably be omitted as we have all the info we need in the ploidy column).
I would like to visualize them on a chromosomal plot. The only tool I've managed to make work with my data is CNANorm (http://www.bioconductor.org/packages/devel/bioc/vignettes/CNAnorm/inst/doc/CNAnorm.pdf), but it insists on analyzing my already analyzed data, so I can't get the results I need. Here's my CNANorm output, although the data points are all wrong because of the extra analysis. ![enter image description here][1]

Any help would be appreciated, please keep in mind that I don't have that much R experience, and that I only need a viewer for the time being.
Thanks everyone :)

[1]: https://dl.dropboxusercontent.com/u/13728657/TEMP/Public%20-%20Don%27t%20overwrite/CNV%20for%20PR.png

Answer: You can plot your data like this. CNVs are plotted as segments with a point corresponding to their midpoint. Chromosomes are separated out into facets.

library("ggplot2")

chr <- c("chr1", "chr2", "chr2")
start <- c(86000000, 70250000, 203050000)
end <- c(117150000, 70500000, 204650000)
center <- start + (end - start)/2
ploidy <- c(1, 3, 3)

ploidy_df <- data.frame(chr, start, end, center, ploidy)

ggplot(ploidy_df, aes(x=center, y=ploidy)) +
geom_point() +
geom_segment(aes(x=start, y=ploidy, xend=end, yend=ploidy, colour="segment")) +
geom_hline(y=2, linetype=2) +
facet_wrap(~chr) +
xlab("Position") +
theme_bw()
![enter image description here][1]

[1]: https://i.imgsafe.org/39e3a1b.jpeg
Question: I usually use `bedtools intersect` to find overlapping regions of bed files, but it seems like this tool can only output overlap between a pair of files. I need something that can do the following from many bed files and **only** report regions contained in all of the bed files.

Example input from 4 separate bed files:

```
chr1 50 100
chr1 60 120
chr1 30 90
chr1 50 90
```

Desired output:

```
chr1 60 90
```

Any tools for this? Maybe I should just `cat` all the bed files together and merge them?

Answer: https://www.biostars.org/p/172566/ (`bedops` option)

There is also

bedtools multiinter

Tool: bedtools multiinter (aka multiIntersectBed)
Version: v2.26.0
Summary: Identifies common intervals among multiple
BED/GFF/VCF files.

Usage: bedtools multiinter [OPTIONS] -i FILE1 FILE2 .. FILEn
Requires that each interval file is sorted by chrom/start.
Question: I have a matrix (or data frame, Excel spreadsheet). It is populated with alphanumeric identifiers (ATG numbers if you are familiar with plants; ex. At1G45623). The first column has only 1 occurrence of each identifier (no duplicates). Each row, after the first column, has a variable number of these alphanumeric identifiers. Anywhere from 2 to 350 in each row. A simplified version is below:

z a b c d e f g h i j

y a b c z w h n q i j

w j i h a f d c b g e p

x i j a b c d h e f m k l o u p s t v

So, each row is like an array where the name of the array is the value in the first column.

I would like to cluster these 'arrays' so that row z (array z) and row w cluster together and not row z with row y. Disregarding the first column, rows z and w have the same group of letters (except for 1 letter), but in a very different sequence. Rows z and y have similar sequences, but somewhat different letters. The method I am searching for would cluster rows z and w together before z and y.

The clustering techniques I am familiar with all take the sequence of the values into account. I am looking for a method that would disregard the sequence and just consider the contents of the row.
This is further complicated by the fact that the rows have very variable numbers of values, and that is why I have included row x just to remember that the size of each 'array' is variable. Does anyone know of anything in Perl, Python or R that could help ?

Answer: Compute the Jaccard index (or any other suitable measure of similarity over sets, see the [sets][1] package in R) and use a clustering algorithm.

[1]: https://cran.r-project.org/web/packages/sets/index.html