cannin/biostars_qa huggingface dataset:
Question: I have a narrowPeak file that is in hg38 format and I need to convert it to hg19 to use it in finding differential looping of HiChIP pipeline.
My general idea was to convert the narrowpeak file into bed file by the following command

cut -f 1-6 MCF10A_H3k27ac_hg38.narrowPeak > MCF10A_H3k27ac_hg38_edited.bed

and using this bed file to liftover to hg19 but I will be losing the remaining narrowpeak data. But again I can't add the remaining narrowpeak file columns to my converted bed file as some of the records are failed to convert in the liftover tool.

Is there a better way to solve this problem?

Thank you!

Answer: Hello,

narrowPeak is an extended bed format. In this case the first 6 fields conform to bed, and then there are an additional 4 fields: https://genome.ucsc.edu/FAQ/FAQformat.html#format12

You can pass the `bedPlus=6` argument to the liftOver util in order to preserve all columns in the file.

Some of this you may be familiar with, but you would need the hg38 to hg19 chain file: https://hgdownload.soe.ucsc.edu/goldenPath/hg38/liftOver/hg38ToHg19.over.chain.gz

The liftOver util (corresponding to your operating system): http://hgdownload.soe.ucsc.edu/admin/exe/

Then you can invoke liftOver on the narrowPeak files directly:

cat exampleHg38.narrowPeak
chr1 9356548 9356648 . 0 . 182 5.0945 -1 50
chr1 9358722 9358822 . 0 . 91 4.6052 -1 40
chr1 9361082 9361182 . 0 . 182 9.2103 -1 75

./liftOver -bedPlus=6 exampleHg38.narrowPeak hg38ToHg19.over.chain exampleHg19.narrowPeak unmapped.txt
Reading liftover chains
Mapping coordinates

cat exampleHg19.narrowPeak
chr1	9416607	9416707	.	0	.	182	5.0945	-1	50
chr1	9418781	9418881	.	0	.	91	4.6052	-1	40
chr1	9421141	9421241	.	0	.	182	9.2103	-1	75

If you have any follow up questions, our public help desk can always be reached at genome@soe.ucsc.edu. You may also send questions to genome-www@soe.ucsc.edu if they contain sensitive data. For any Genome Browser questions on Biostars, the UCSC tag is the best way to ensure visibility by the team.
Question: Hi,
I have two files:

File1:
chr1 73270408 73271242 XLOC_000129
chr1 83659841 83660368 XLOC_000153
chr1 163894423 163894714 XLOC_000257
chr1 55485657 55486215 XLOC_000374
chr1 155285623 155287724 XLOC_000553
chr1 155288475 155288592 XLOC_000553
chr1 155289366 155290533 XLOC_000553
File2:
chr1 53933832 53934533 TCP11L2
chr1 53935094 53935266 TCP11L2
chr1 53935961 53936136 TCP11L2
chr1 53937602 53937789 TCP11L2
chr1 53938748 53938884 TCP11L2
chr1 53939440 53939660 TCP11L2
chr1 53940084 53940204 TCP11L2
chr1 53941131 53941263 TCP11L2
chr1 53942563 53942732 TCP11L2
chr1 53945014 53945137 TCP11L2
chr1 121715553 121715945 GRPR
chr1 121719096 121719447 GRPR
chr1 121736572 121736990 GRPR
How can I want to extract genes from File2 that are located 10000 bp upstream and 10000 bp downstream of XLOC_* in File1?

Answer: Via [BEDOPS][1], you can use `bedmap` to get an answer, depending on whether you want gene names and intervals, or just gene names.

First, sort your files with BEDOPS `sort-bed`:

$ sort-bed File1.unsorted.bed > File1.bed
$ sort-bed File2.unsorted.bed > File2.bed

You only need to do this sort step once.

If you want *gene names and intervals* from File2, overlapping elements from File1 that are padded by 10k up- and downstream:

$ bedmap --range 10000 --echo-map --multidelim '\n' --skip-unmapped File1.bed File2.bed | sort-bed - | uniq - > answer.bed

If you just want *gene names*:

$ bedmap --range 10000 --echo-map-id File1.bed File2.bed | sort - | uniq - > answer.txt

If you want to know which elements from File2 overlap 10k-padded elements from File1:

$ bedmap --range 10000 --echo --echo-map File1.bed File2.bed > answer.bed

In this case, each line of `answer.bed` is an original element from File1, along with elements from File2 that overlap the 10k-padded File1 element.

See the [documentation][2] for more help, or run `bedmap --help`.

[1]: https://github.com/bedops/bedops/
[2]: http://bedops.readthedocs.io/en/latest
Question: The plot is from this publication, see page 2 the bottom of the plot with ACTG letters in different colours and sizes.

Could anyone clarify in layman's terms what it represents?

More importantly, how do we plot it using R or any other tools?

Here is the plot, if you can't access the paper:

Found related post at SO: Plotting a “sequence logo” using ggplot2?

Answer: That type of plot is called Sequence Logo. The X axis represent the nucleotide positions of a sequence, and the Y the entropy of each possible nucleotide. The size of each letter indicates how frequently it is found at that position: for example position 8 is almost always an A, while position 5 can be a T or C.

You can create sequence logos from a web interface using a popular application called WebLogo, or using R packages such as SequenceLogo.
Question: Can anybody help me to find the dataset for s. pombe on BioMart? And also some help on how to use `makeTranscriptDbFromBiomart` to create TranscriptDB?

cheers,

Answer: Looks like you figured out another way of getting what you needed, but, for the record, here is the answer to your question:

S pombe is at http://fungi.ensembl.org/index.html

The biomart is here: http://fungi.ensembl.org/biomart/martview/248a3d2deec76fa7be1e94e32b3972df

Access it using BioConductor's GenomicFeatures as follows. Note the warnings....

library(GenomicFeatures)
library(biomaRt)

txdb<-makeTranscriptDbFromBiomart(
,biomart ="fungi_mart_22"
,dataset = "spombe_eg_gene"
,host="fungi.ensembl.org"
)

Download and preprocess the 'transcripts' data frame ... OK
Download and preprocess the 'splicings' data frame ... OK
Download and preprocess the 'genes' data frame ... OK
Prepare the 'metadata' data frame ... OK
Make the TranscriptDb object ... OK
Warning messages:
1: In .normargSplicings(splicings, transcripts_tx_id) :
no CDS information for this TranscriptDb object
2: In .normargChrominfo(chrominfo, transcripts$tx_chrom, splicings$exon_chrom) :
chromosome lengths and circularity flags are not available for this TranscriptDb object

> transcriptsBy(txdb)
GRangesList of length 7017:
$SPAC1002.01
GRanges with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
|
[1] I [1798347, 1799015] + | 510 SPAC1002.01.1

$SPAC1002.02
GRanges with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
[1] I [1799061, 1800053] + | 511 SPAC1002.02.1

$SPAC1002.03c
GRanges with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
[1] I [1799915, 1803141] - | 2075 SPAC1002.03c.1

...
<7014 more elements>
---
seqlengths:
I II III MT MTR AB325691
NA NA NA NA NA NA
Question: Hello,

I'm trying to automate clustalw workflow, but as result I need 3 output files .fasta .nexus and .dnd. I've tried many parameters but as result i get always only 2 files .fasta or .nexus and .dnd. The command I'm trying to enter into linux terminal :

clustalw -infile=third.fna -align -type=DNA -score=PERCENT \
-outfile=out.fasta -output=FASTA -outfile=out.nexus \
-output=NEXUS -maxdiv=45

Do anyone have idea how should I fix my command either it's just impossible for clustalw to ouput that amount of files at once?

Answer: I've figured it out. There's no way to produce extra output file from ClustalW like [jrj.healey][1] said. More over there are two simply software's that allow us to convert .fasta file to .nexus format.

pip install seqmagick
seqmagick convert --output-format nexus --alphabet dna inFile.fasta outFile.nexus

also suggestion by [a.zielezinski][2] with:

sudo apt-get install readseq
readseq inFile.fasta -all -format=NEXUS -output=outFile.nexus

Thanks for help!

[1]: https://www.biostars.org/u/20598/
[2]: https://www.biostars.org/u/4700/
Question: I have vcf data from 95 people - 63 Cases and 32 controls. I filtered out possible pathogenic variations and now i want to know if variation in Cases samples apper significantly more often then in controls parth of sample.

I decided for Fst statistics done with vcftools. However I am not still 100% what am i doing.

Is there betters test to test for differenciations of two population? Is there any better tools then vcftools?

Answer: you could use `bcftools contrast` https://samtools.github.io/bcftools/howtos/plugins.html

> runs a basic association test, per-site or in a region, and checks for novel alleles and genotypes in two
>
> groups of samples. Adds the following INFO annotations:
>
> PASSOC .. Fisher’s exact test probability of genotypic association (REF vs non-REF allele)
Question: Dear all,

I have a reference that looks like this
```
AAAAATTTTTTTNNNNNNGGGGGGCCCCCC
```
and on this reference reads have been mapped

```
ref: AAAAATTTTTTT-**NNNNNN-**GGGGGGCCCCCC

read1: AAAAATTTTTTT-**ATGCAA**-GGGGGGCCCCCC
....
read2: AAAAATTTTTTT-**CCGTAA**-GGGGGGCCCCCC
```

and then I extracted the bam file.
I want to extract the sequences that have been mapped the reference in the positions 13:18.
For this I have used the following command
```
samtools view -b -h file.sorted.bam Reference_barcodes:13-18 > test.bam
```
Now I want from this bam to isolate in a fasta file the exact sequences that have mapped the 13-18 and my data to look
like this

```
ref: NNNNNN
read1:ATGCAA
read2:CCGTAA
```

Answer: Finally only this gives the answer I want:

```
library(GenomicAlignments)
stack <- stackStringsFromBam("test_4_seqs_sorted.bam", param=GRanges("Reference_barcodes:54-78"))
```
Question: I have a couple of questions about comments or documenting in CWL.

1a. What is the difference between 'doc' and 'label' in CWL?

1b. What are the special cases in which one is better than the other?

2a. Are line that start with hash tags, also comments?

2b. if no, what is the use of hash tags? if yes, whats the difference between a hash tag, 'doc' and 'label'?

Answer: Hello kevin.o.oluoch and thank you for your question.

1a:

From http://www.commonwl.org/v1.0/CommandLineTool.html#CommandLineTool we learn that

> `label` […] A short, human-readable label of this process object.

and

> `doc` […] A long, human-readable description of this process object.

1b. `doc` is really for longer documentation, possible several paragraphs. `label` shouldn't be more than a sentence (if not just a fragment)

2a. hash tags (YAML comments) are comments visible only to viewers of the raw CWL description. They don't get represented in the data model.

Example: https://github.com/ProteinsWebTeam/ebi-metagenomics-cwl/blob/ca6ca61/tools/FragGeneScan1_20.cwl#L35

inputs:
sequence:
type: File
format: edam:format_1929 # FASTA

2b. I personally use the in line comments to leave notes for myself and others who are editing or reading the raw CWL description.

I hope you found this explanation helpful!
Question: Hi guys, sorry if I look dumb (I do feel that way now).

I am going through info pages, citations and database metrics of Blast since this morning because I need to find these infos:

- Number of sequences in the non-redundant protein database (nr)
- Number of sequences in the non-redundant nucleotide database (nt)

I am a little desperate. I don't think that the number of sequences in NCBI-Genbank and the number of sequences in NCBI-Protein are the same as in **nr** and **nt**, because they're probably redundant. Am I wrong? Do you know where could I get such numbers?

**EDIT:** I've been suggested a **workaround** elsewhere so I post it for future readers.
The database folder contains the information in a file.

- nt.nal for the nucleotide nr
- nr.pal for the protein nr

Answer: nt.nal and nr.pal files contains info about the number of sequences (NSEQ).
Details from these files for our local BLAST database version:

#
# Alias file created 06/24/2017 23:11:41
#
TITLE Nucleotide collection (nt)
DBLIST "nt.00" "nt.01" "nt.02" "nt.03" "nt.04" "nt.05" "nt.06" "nt.07" "nt.08" "nt.09" "nt.10" "nt.11" "nt.12" "nt.13" "nt.14" "nt.15" "nt.16" "nt.17" "nt.18" "nt.19" "nt.20" "nt.21" "nt.22" "nt.23" "nt.24" "nt.25" "nt.26" "nt.27" "nt.28" "nt.29" "nt.30" "nt.31" "nt.32" "nt.33" "nt.34" "nt.35" "nt.36" "nt.37" "nt.38" "nt.39" "nt.40" "nt.41" "nt.42" "nt.43" "nt.44" "nt.45" "nt.46" "nt.47"
NSEQ 43107468
LENGTH 148603982942
#
# Alias file created 06/09/2017 09:55:12
#
TITLE All non-redundant GenBank CDS translations+PDB+SwissProt+PIR+PRF excluding environmental samples from WGS projects
DBLIST "nr.00" "nr.01" "nr.02" "nr.03" "nr.04" "nr.05" "nr.06" "nr.07" "nr.08" "nr.09" "nr.10" "nr.11" "nr.12" "nr.13" "nr.14" "nr.15" "nr.16" "nr.17" "nr.18" "nr.19" "nr.20" "nr.21" "nr.22" "nr.23" "nr.24" "nr.25" "nr.26" "nr.27" "nr.28" "nr.29" "nr.30" "nr.31" "nr.32" "nr.33" "nr.34" "nr.35" "nr.36" "nr.37" "nr.38" "nr.39" "nr.40" "nr.41" "nr.42" "nr.43" "nr.44" "nr.45" "nr.46" "nr.47" "nr.48" "nr.49" "nr.50" "nr.51" "nr.52" "nr.53" "nr.54" "nr.55" "nr.56" "nr.57" "nr.58" "nr.59" "nr.60" "nr.61" "nr.62" "nr.63" "nr.64" "nr.65" "nr.66" "nr.67" "nr.68"
NSEQ 124946261
LENGTH 45844202229
Question: Dear all,

I tried to use Picard MarkDuplicates to mark duplicates. I found `METRICS_FILE` need to be set. However, I know nothing about the "metrics"about my data. All I know is its mapping rate, etc. What should I look into and how should I set the matrics file please? Thank you.

Yingzi

Answer: Hi, my syntax looks like:

picard -Xmx"$memory"g -Djava.io.tmpdir=$PWD/tmp -XX:ParallelGCThreads=$cpu MarkDuplicates I= ${p%.bam}.bam O= ${p%.bam}_rmdp.bam METRICS_FILE= ${p%.bam}_rmdp.txt CREATE_INDEX=true USE_JDK_DEFLATER=true USE_JDK_INFLATER=true

Metrics should give you an statistics about PCR duplicates.
Question: First off, I apologize for posting this "old" inquiry. I know similar inquiries were put out before, but I am hoping to find new solutions to this inquiry.

I am trying to find the positions of a list of SNPs (given the rs#). Need a "new" solution.

What I have tried so far -

(a) sending a batch query to dbSNP at NCBI, which worked well in the past, but today ~10 hr after sending the batch query, no return of result yet ( is the server down ??)

(b) downloading all dbSNP positions from Biomart and hoping to do some "intersection" to find the positions for specific rs#. The download somehow was terminated prematurely (first download took ~ 1+ hr).

(c) downloading cruzdb. cruzdb was suggested as a solution in one of the earlier posts. I read the document and still could not run it - my apology ! (does running cruzdb require an understanding of python? which I currently don't possess) Having to say it though, in contrast to cruzdb doc, I had better luck with vcftools and plink thanks to their "more friendly" documents.

Is there any other solutions that allow non-bioinformaticians to find answers to this task (i.e. positions for a list of SNPs)?

I certainly hope to get some useful responses, but It's understandable if the admin chooses to close this thread (due possibly to "duplication of questions"). Thank you

Answer: Hi, as suggested by Ashutosh Pandey, you can exploit [UCSC table browser][1]. Select genome and assembly of interest and from **group** menu select *Variation*. From **track** menu select *All SNPs(142)*. Paste or upload your rs ids using buttons at **identifiers (names/accessions)**. To export your results, select *selected fields from primary and related tables* from **output format** then click get output. In the next step you can select fields of interest (*i.e.* input rds, chromosome, genomic position) that will be included in final output table, click *get output* to retrieve it.

[1]: https://genome.ucsc.edu/cgi-bin/hgTables
Question: Hi,
I have a set of piRNAs in fasta format and I want to find their potential targets.
According to literature nucleotides in position 1-10 in the piRNAs sequences should constitute the seed and do not tolerate any mismatch while nt 11-20 tolerate up to 3 mismatches.
Bowtie can specify the length of the seed, but I'm having troubles finding the right combination of parameters to translate the piRNA target pairing logic into a single bowtie command line.

Is there any solution to get this done with a single alignment run or do I have to run a first alignment guided by the seed 1-10nt, and successively another alignment for the 11-20 region allowing for up to 3 mismatches?

Thank you in advance for your help!

Answer: How about running the alignment tolerating 3 mismatches anywhere, and then filtering the results using the CIGAR string ? (You would need an aligner that produces CIGAR strings distinguishing matches (`=`) and mismatches (`X`); I do not remember if `bowtie` does).
Question: I am trying to use samtools sort function to sort a bam file (named as PT2) (size 180 megabite), which is not large. The sorting never end and no error message appeared. I have checked the memory of the VM Ubunto using (meminfo file in the proc folder on the home account). Also I have around 80 gig disc space

First code (extremely slow) :

samtools sort -o outputs/PT2_sorted.bam outputs/PT2.bam

second code (very slow):

samtools sort -o outputs/PT2_sorted.bam outputs/PT2.bam

What do you think is a reason behind this slow working ? any advice ?

Thanks

Answer: That was resolved. The point was to increase the VM memory from 1024 to around 1800. This enabled running the command (even without adding `-m` or `-T`)

Thanks
Question: Hi Guys,
I am new to the forum to come directly with a question / problem.

I want to write a program that counts the patterns like this.
I have given such sequences as an example:

> - A: AAAAAAAAA
- B: AAAAACAAG
- C: GAAAACGAA
- D: AAAAAAAAT
>

A patten is a group of four letters of the four groups. As an example, the first Patten would be "AAGA". Now I wonder how many times that happens in my entire Sequenze. So far I have only found programs that are looking for a given pattern. Does anyone have an idea how to program this or do you already know where this was done?

I'm very thankful for your help!

Best regards

Answer: Something like this ?

#### Input fasta file : test.fa

>A
AAAAAAAAA
>B
AAAAACAAG
>C
GAAAACGAA
>D
AAAAAAAAT

#### Python counter

from Bio import SeqIO
import pandas as pd
from collections import Counter

df = pd.DataFrame()
for record in SeqIO.parse("test.fa", "fasta"):
df[record.id] = list(str(record.seq))
df = df.transpose()

It returns a dataframe like this :

0	1	2	3	4	5	6	7	8
A	A	A	A	A	A	A	A	A	A
B	A	A	A	A	A	C	A	A	G
C	G	A	A	A	A	C	G	A	A
D	A	A	A	A	A	A	A	A	T

Then you can count column string occurence with Counter

c = Counter([df[i].str.cat() for i in df])
# Output : Counter({'AAGA': 2, 'AAAA': 5, 'ACCA': 1, 'AGAT': 1})

Note :
It only works if you have few sequence. Loading the dataframe with a huge fasta will fail with low memory.
You also may want to see [Position Weight Matrix][1]

[1]: https://en.wikipedia.org/wiki/Position_weight_matrix
Question: I'm using biopython to perform an Entrez esearch of the NCBI taxonomy database. The term is the exact string as it appears in the database.

Entrez.read(Entrez.esearch(db='Taxonomy', term="Cupriavidus sp. HPC(L)"))

returns

{'QueryTranslation': 'Cupriavidus sp. HPC[All Names] AND (L[All Names])', 'IdList': [], 'Count': '0', 'RetStart': '0', 'ErrorList': {'FieldNotFound': [], 'PhraseNotFound': ['Cupriavidus', 'sp.', 'HPC', 'L']}, 'RetMax': '0', 'TranslationSet': [], 'WarningList': {'PhraseIgnored': [], 'QuotedPhraseNotFound': [], 'OutputMessage': ['No items found.']}}

This should return information on TaxID 1217418.

Thanks in advance.

Answer: SIMPLE AND OBVIOUS SOLUTION:

Enclose the search term in double quotes to make it literal.

Entrez.read(Entrez.esearch(db='Taxonomy', term="\"Cupriavidus sp. HPC(L)\""))

correctly returns

{'RetStart': '0', 'TranslationStack': [{'Count': '1', 'Field': 'All Names', 'Term': '"Cupriavidus sp. HPC(L)"[All Names]', 'Explode': 'N'}, 'GROUP'], 'IdList': ['1217418'], 'TranslationSet': [], 'Count': '1', 'RetMax': '1', 'QueryTranslation': '"Cupriavidus sp. HPC(L)"[All Names]'}
Question: Hi,

I'd like to automatically download the tar files from the following links, does anyone knows how to do that?

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39977

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36375

Kind of like using wget...

Thanks!

Answer: Never mind, just figured it out. I can go to the ftp site and do wget:

ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE39nnn/GSE39977/suppl/

wget 'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE39nnn/GSE39977/suppl/GSE39977_RAW.tar'

Thanks!
Question: Hi,

I am trying to write a r code to calculate the pearson correlation for the first row vs the rest of the rows. (1vs2, 1vs3,1vs4,1vs5 so on). I am not sure if apply,lapply or sapply will work for this. Also when I try to loop it says x and y must be of same length, but here always my x is going to row 1 and y is 2:19725.

How can make it workable in R

Thanks,

Answer: R function `cor()` automatically calculates the correlation between the columns of matrix, so you can transpose your matrix and calculate it like this:

```
# Creates a random matrix with 5 rows and 6 columns
mat=matrix(rnorm(30),5,6)
# Computes the correlation between the first row and the other 4
cor(mat[1,],t(mat[2:5,]))
```
Question: I have seen some other posts on here but after going through them Im not finding a solution. I am using limma to do DE analysis on 763 patient microarrays, 4 groups total. My expression object is a matrix with rownames as genes and colnames as patient ID.

When I try to run the code I get the error:

> fit <- lmFit(as.numeric(Merge_DF_Avg2), design)
Error in lmFit(as.numeric(Merge_DF_Avg2), design) :
row dimension of design doesn't match column dimension of data object

The dimensions seem correct:

> dim(design)
[1] 763 4

> dim(Merge_DF_Avg2)
[1] 20341 763

> class(Merge_DF_Avg2)
[1] "matrix" "array"

Here is the code:

#Data wrangle for limma
>Merge_DF_Avg <-readRDS("C:/Users/12298/Desktop/Data_Analytics/Taylor_2017/Finalized_Expression_Array_avgs.rds")
>Merge_DF_Avg2 <- Merge_DF_Avg %>% na.omit() %>% pivot_wider(values_from=avg, names_from=Gene_Name) %>% t() %>% janitor::row_to_names(row_number = 1)
>Patient_Cat <- as.vector(as.numeric(Merge_DF_Avg2[1,]))
>Merge_DF_Avg2 <- Merge_DF_Avg2[-c(1,2),] #Taking out patient cat and unidentified gene rows

#limma design
>design <- model.matrix(~ 0 + factor(Patient_Cat))
>colnames(design) <- c("SHH", "Group3", "Group4", "WNT")
>fit <- lmFit(as.numeric(Merge_DF_Avg2), design)

I am using the example code provided by the limma package pdf to guide me.

Answer: I think `as.numeric(Merge_DF_Avg2)` in `lmFit` converts a character matrix to a vector. To convert the numeric vector back to matrix you could do (check it's ok!):

matrix(as.numeric(as.matrix(Merge_DF_Avg2)), ncol= ncol(Merge_DF_Avg2))