cannin/biostars_qa huggingface dataset:
Question: Hi,

I am new to Bioinformatics.I am trying to use GATK tool for finding SNP and indels but the problem is that the documentation seems to be complex for a beginner like me with so many tools to start with including IndelRealigner,MuTect,Haplotypecaller..The more I read the more confused I am getting.I also tried using Indel religner which is the first step in the process of identification of indels(correct me if wrong) but I got an error saying that you don't have the intervals file in the targetIntervals in the following command.

java -jar GenomeAnalysisTK.jar \
-T IndelRealigner \
-R reference.fasta \
-I input.bam \
-known indels.vcf \
-targetIntervals intervalListFromRTC.intervals \
-o realignedBam.bam

The error message is:Could not read file /home/aditya/Bioinfotools/gatk/intervalListFromRTC.intervals because The interval file does not exist.
I don't know where to get the intervals file from.

But please guide me overall in understanding GATk.

Answer: Realignment is a two-step process using RealignerTargetCreator and then IndelRealigner. See detailed instructions here:
https://www.broadinstitute.org/gatk/guide/article?id=38

The rest of the best practices tutorial will be helpful too.
Question: I've recently upgraded to IGB v8.3.0 on OSX, but I don't seem to have the option for choosing different RAM version as there used to be.

Anyone know how to change the RAM requirement? Currently only able to access 1GB.

Answer: Nowlan has already done a great job explaining the new default behavior, but I will add a note for power users. There is a "vmoptions" file in your IGB installation directory which contains a commented out jvm argument for setting the max heap size. If you want to manually control the heap size you simply need to uncomment this line and adjust it to your needs.

- On Mac
- `/Applications/Integrated\ Genome\ Browser/IntegratedGenomeBrowser.app/Contents/vmoptions.txt`
- On Linux/Windows
- `${IGB_INSTALLATION_DIR}/IntegratedGenomeBrowser.vmoptions`
Question: I have a list of HLA alleles, but I need to get these IDs in rsID format, or at least the chromosome and position.

I've tried searching on SNP-DB and SNPedia, and old-fashioned internet search, but nothing is coming up.

PharmGKB doesn't show the necessary position or rsID.

How can I connect an HLA allele, for example HLA‐C*07:29 as an rsID or chromosome and position?

Answer: HLA are gene haplotypes as in defined by the thousands of nucleotides that constitute an HLA gene. I could see why they would not have rsids and positions like single nucleotide variants and short insertion deletions. Happy to be shown otherwise though I didnt research this at all
Question: I know for the reverse procedure, there are lots of methods available. But how can I change my fasta file with sequences like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRQQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRVKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

to like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
QQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
VKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

I know its just about inserting line breaks in each sequence, but is there any way that I can check whether the fasta file is single line fasta or not, and if is, add line breaks in each sequence..

Ps. I am using windows..

Thanks for your consideration

Answer: [seqmagick](https://github.com/GuangchuangYu/seqmagick):

```
library(seqmagick)
fa_read(INFILE) %>% fa_write(OUTFILE, type = "interleaved")
```
Question: HI,

Do I have to remove repeats in the newly sequenced genomes ?
The genomes will be used for Phylogenomics study.

Suggestions please!

Thanks

Answer: yes and no.

it is generally always good to mask the repeats in (new) genomes. Don't remove them though! but rather soft (or hard) mask them.

On the other hand, those repeats should not interfere much with phylogenomics analyses, as those kind of studies are primarily based on 'normal' protein coding genes. In that sense it is important to make sure the repeat(proteins) do not end up in your proteome for a certain organism.
Question: Hi,

I am new to variant analysis. I have mutation induced RNAseq samples data which has been aligned to reference transcriptome and obtained a bam file which has been sorted. For every transcript I would like to count the number of mutated reads aligned to reference transcript. For an example, number of reads on which base C has been mutated to A at every position along the transcript.

Is it possible to achieve this using samtools and vcf tools.

Kindly guide me

Answer: [bam-readcount][1]?

[1]: https://github.com/genome/bam-readcount
Question: Dear All,

I have to deal some small molecules screening for my protein-drug docking studies. Please suggest me few small molecule databases ..also the free software can use to screen them...thank you

Answer: As you suggested, you will find a smaller but usefully target mapped compound and drug set at http://www.guidetopharmacology.org/

You should probably start docking with PDB structures anyway so you can get a small set from the above by using the PubChem query

https://www.ncbi.nlm.nih.gov/pccompound?term=%28%231%29%20AND%20pccompound_structure[filt]%20AND%20%22IUPHAR/BPS%20Guide%20to%20PHARMACOLOGY%22[SourceName]&cmd=DetailsSearch=872
Question: I used bowtie2 to map my reads to a reference genome and get the SAM files returned. I am wondering if there is a way to get the length of the identical mapped nucleotide length. Basically I want to get a hamming distance or identical mapping vs total read length ratio from the SAM files.

Answer: find out the answer in another post: https://www.biostars.org/p/7331/
Question: Hey,

I have to do some variant calling in a plant organism (diploid)

I have selected the programs I want to test before I use my own data. (Freebayes, Samtools, SOAPsnp and Cortex_VAR.)
However i am struggling to come up with a reference to compare the tools against.
First i want to test the software before i can use it on my own data.

Is there a benchmark i can use or is it better to simulate my own reference? and use the tools accordantly

Regards,
Covux

Answer: use 'Genome in a Bottle' data ? http://jimb.stanford.edu/giab/

e.g: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/
Question: I'm using Pfam HMMs to identify domains in my *T. cruzi* proteins with the HMMER program.

I kinda understand that a HMM is build by a family multiple sequence alignment, by calculating probabilities of states.

But how the hmmscan works to classify a new protein by families HMMs? Where can I read more about this?
Academic references would be appreciated.

Answer: In this context, HMMs are numerical representations of sequence alignments that are more information-rich because both residue substitution patterns and gap penalties are treated in probabilistic fashion. Diverse sequences are given higher weight in HMMs, and they also include background residue probabilities in a way that is inversely proportional to the alignment depth.

If you don't have access to Krogh, Durbin & Eddy book, the PDF version is available [**here**][1]. Any of the early Krogh, Eddy and Karplus papers would be good as starting points.

https://www.ncbi.nlm.nih.gov/pubmed/8744772

https://www.ncbi.nlm.nih.gov/pubmed/9918945

https://www.ncbi.nlm.nih.gov/pubmed/9927713

https://www.ncbi.nlm.nih.gov/pubmed/18075166

[1]: https://www.semanticscholar.org/paper/Biological-Sequence-Analysis%3A-Probabilistic-Models-Durbin-Eddy/571f5bbecd3a083a2bb6844f59a3f8cea237252e
Question: **How do you classify the miRNAs? What is the best way to classify miRNAs hypothetical using Blastn-short or RFAM? [doubts]**

Hello everyone. Let's me explain my doubts. I mean about 200 sequences which was predicted by mir-Bag (pre-microRNA) before. Then I used this pre-micros with the miRdup in order to obtain the microRNAs. Now, I want to know which have been listed above. Firstly, I did a `blastn -task blastn-short -evalue 10` (v. 2.2.30) against the mature miRNAs from mirBase.org. I wonder if you use this strategy and what options you usually use like:

```
-evalue 10
-penalty -4?
-word_size 4 or 7?
-ungapped?
-reward 1?
```

Because it is not neccesary to have a full coverage of the query, therefore there are some isomeRs with slightly mismatches. But if the seed region (2nt-7nt) is covered fully then it is the same family miRNA. Isn't it? Or I have to do the RFAM. Regards.

I attach one example and some options of blastn (v. 2.2.30.)

```
query= ENTRY327687_trypsinogen_354_553

Length=22
Score E
Sequences producing significant alignments: (Bits) Value

dme-miR-4941-5p MIMAT0020148 Drosophila melanogaster miR-4941-5p 24.3 0.21
dvi-miR-9538-3p MIMAT0035619 Drosophila virilis miR-9538-3p 22.3 0.83
mmu-miR-7677-5p MIMAT0029868 Mus musculus miR-7677-5p 22.3 0.83
mmu-miR-6922-3p MIMAT0027745 Mus musculus miR-6922-3p 22.3 0.83
gga-miR-6708-5p MIMAT0025818 Gallus gallus miR-6708-5p 22.3 0.83
mtr-miR2590j MIMAT0021334 Medicago truncatula miR2590j 20.3 3.3
mtr-miR2590i MIMAT0021333 Medicago truncatula miR2590i 20.3 3.3
mtr-miR2590h MIMAT0021332 Medicago truncatula miR2590h 20.3 3.3
tca-miR-3813-5p MIMAT0018648 Tribolium castaneum miR-3813-5p 20.3 3.3
mtr-miR2598 MIMAT0013301 Medicago truncatula miR2598 20.3 3.3
oan-miR-181c-3p MIMAT0007059 Ornithorhynchus anatinus miR-181c-3p 20.3 3.3

Query_29 2 TGGCGGTGAGCAGAATAATTG 22
20244 19 GGCGGTGAGCAG 8
34141 12 CAGAATAATTG 22
28818 4 GGTGAGCAGAA 14
26045 13 GGTGAGCAGAA 3
25597 11 GCGGTGAGCAG 21
21071 12 GCAGAATAAT 21
21070 12 GCAGAATAAT 21
21069 12 GCAGAATAAT 21
18383 10 GAGCAGAATA 1
13217 19 GCAGAATAAT 10
7532 10 TGGCGGTGAG 1

Lambda K H
1.37 0.711 1.31

Gapped
Lambda K H
1.37 0.711 1.31

Effective search space used: 4258254
```

Answer: My two cents.

If you want to obtain homologous, I tried once with BLAT modifying the parameter stepSize=2 and tileSize=8 to get all possible alignments. I would do two strategies:

1. Map my precursors to the miRBase precursors to obtain homologous precursors, maybe anything > 80%
2. Map my miRNAs to the miRBase precursors and check if they are inside the corresponding mature sequences allowing an error of 3 nt at the borders. That means, if the mature of that precursor start at the position 23, I will check any of my sequences starting from position 20 to 26. The same for the 3' end.

In both of your case you need to clean your results to get the most score hits, and check them to see if they make sense.
Question: Hi, everyone!

k-tuple frequency can be used in clustering sequences.

I wonder if this method can be used when target sequences are very short like reads from next generation sequencing? Will it be unstable? (I have read some papers and all of them just use k-tuple frequency to catalog different meta-data, but not clustering reads in those meta-data.)

Any hints or paper about this would be helpful! Thank you!

Answer: Yes, they certainly can. A good example would be a newly released piece of software called [Sailfish (Published recently in Nature Biotechnology](http://www.nature.com/nbt/journal/v32/n5/full/nbt.2862.html)) which uses k-mers to do alignment-free isoform assignment with RNA-Seq data.
Question: I am new to bioinformatics and currently learning how to use Bowtie 2. As written in the manual:

> A pair that aligns with **the expected relative mate orientation** and with the expected range of distances between mates is said to align "concordantly". If both mates have unique alignments, but the alignments do not match paired-end expectations (i.e. the mates aren't **in the expected relative orientation**, or aren't within the expected distance range, or both), the pair is said to align "discordantly".

I have read about the basics of paired-end sequencing and orientation (in the molecular biology sense). In summary, my understanding is that in paired-end sequencing we sequence both ends of a DNA fragment at the 5' end and the 3' end (we call them mate 1 and mate 2) and by knowing the expected length between the mates we can better align the fragment to a reference genome.

My question is, what is it meant by two mates having an expected relative orientation? If am using Bowtie 2 and giving it a file with all the first mates and another one with all the corresponding second mates, how can two mates align without first having the expected orientation?

Answer: In Illumina paired-end sequencing the two two mates will be sequenced in the opposite direction to one another (see my ASCII art below), so what you expect from this is that the reads should be pointing into each other, which would be the expected orientation. The reads can be in other orientations if:

1. say you were mapping to a reference from a related organism which had an inversion in its genome relative to what you sequenced;
2. The reads mapped to a palindromic region and the read mapper could not determine the correct orientation
3. There is a problem with the assembly that you are mapping to (for example errors in a de novo assembly)

-------> Read 1
========================= DNA fragment
<--------- Read 2

Also note that other sequencing technologies produce reads in different orientations. For example, Illumina mate-pair creates reads that point away from each other rather than into each other. I don't know how bowtie handles these cases where other orientations are expected
Question: Hi all,

I'm trying to write some code to generate mean read length data from a range of fastq files.
awk '{if(NR%4==2) print NR"\t"$0"\t"length($0)}' HG1.fastq > readLength.txt

i've got as far as here from looking through other posts and trying to improve but i'm stuck on a couple of things. This command only works on a single file and will report the length of each read within that file separately.

I want to run a single command so the mean and Standard Dev of read lengths from all .fastq files within a folder are reported in a single .txt file, one sample per line. I gues SD might be difficult to calculate in a command so even just the mean read length.

e.g.the first 5 files in my folder are:
ru1.fastq
ru2.fastq
hg3.fastq
hg25.fastq
ru7.fastq

obviously i'm a bit of a novice at this so all help would be appreciated !!

thanks a lot

Answer: Using awk:

awk 'BEGIN { t=0.0;sq=0.0; n=0;} ;NR%4==2 {n++;L=length($0);t+=L;sq+=L*L;}END{m=t/n;printf("total %d avg=%f stddev=%f\n",n,m,sq/n-m*m);}' *.fastq
Question: Dear all,

I have some questions regarding demultiplexing Illumina BCL files. First question is if Illumina´s Undetermined Fastq files contains barcodes?

Second question - If I would like to provide demultiplexing and keep my barcodes in fastq files (no in the header like CASAVA provide, but inside the read sequence) - does anybody have any experience how to do that?

Thak you for any sharing your experience.

Answer: Yes, the undetermined files have the barcodes in the header of each read.

I haven't a clue how you would include read #2 into read #1, though perhaps playing with --use-bases-mask would do the trick (presumably by adding the index length to that of read1). It would make more sense to me to just write read 2 as separate file (so paired end experiments will have all three reads rather than just two).
Question: Dear all,

please would you let me know if there are any Hi-C databases that would have collected the 3C interactions reported by individual laboratories, or big consortia such as ENCODE or 4Dnucleome ? thanks a lot,

-- bogdan

Answer: [This list](https://github.com/mdozmorov/HiC_data) on github is a pretty good one in terms of datasets/resources. [3DIV](http://kobic.kr/3div/) might be a little more what you're after though.
Question: I used this R script of rnaseqdata to generate a Volcano plot:

but there are 3 outliers point deforming the graphic shape. How can I remove it ?

```
resultado_table=rbind(data.frame(results(dds, contrast=c("condition", "Controle" , "Infectado"))))
resultado_table$Variable[1:39179] <- "Controle_vs_Colombiana"

##Highlight genes that have an absolute fold change > 2 and a p-value < Bonferroni cut-off
resultado_table$threshold = as.factor(abs(resultado_table$log2FoldChange) > 1 & resultado_table$padj < 0.05)
write.table(resultado_table, file="./GENES_DE")

volcano_plot <- ggplot(data=resultado_table, aes(x=log2FoldChange, y=-log10(padj) , colour=threshold)) +
geom_point(alpha=0.4, size=1.75) +
xlab("\n log2 fold change") +
ylab("-log10 p-value adjusted \n") +
theme_bw() +
theme(axis.title.y = element_text(face="bold", size=16),
axis.title.x = element_text(face="bold", size=16),
axis.text = element_text(size=12),
legend.title =element_blank() ,
legend.text = element_text(size = 12)) +
facet_wrap(~ Variable, ncol=3) +
theme(strip.text = element_text(size=12, face="bold"))
volcano_plot
```

Answer: If you really want to remove the outliers you could remove the values for the x and/or y axis the data points outside a given quantile range. E.g. to remove values less than 0.01 and greater than 0.99 quantile:

```
qqx<- quantile(resultado_table$log2FoldChange, p= c(0.01, 0.99))
qqy<- quantile(resultado_table$padj, p= c(0.01, 0.99))

resultado_table2<- resultado_table[
(resultado_table$log2FoldChange > qqx[1] & resultado_table$log2FoldChange < qqx[2]) &
(resultado_table$padj > qqy[1] & resultado_table$padj < qqy[2]), ]
```

Now use `resultado_table2` for plotting.
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: You could use a [Python regular expression][1] with a [lookahead][2] for your [kmer][3] of interest:

>>> import re
>>> a = 'AAAAACAAGAAA'
>>> [x.start() for x in re.finditer('(?=AAA)', a)]
[0, 1, 2, 9]
>>> [x.start() for x in re.finditer('(?=AAZ)', a)]
[]

The `len()` property of the list of "hit" start positions is the number or frequency of the kmer in question.

Consider also what you might need to do in Python, if you need to search the [reverse complement][4] of the query sequence.

Another possible approach would be to write a [Python generator][5] of all contiguous substrings of a query sequence (and, possibly, its reverse complement). Using a [generator][6] would reduce the memory overhead of making a list from all substrings.

One can iterate or walk through the output of a generator, one element at a time. You might compare the generator output with your kmer of interest and increment a counter, if there's a match.

[1]: https://docs.python.org/2/library/re.html
[2]: https://www.regular-expressions.info/lookaround.html
[3]: https://en.wikipedia.org/wiki/K-mer
[4]: https://bioinformatics.stackexchange.com/a/3595/776
[5]: https://stackoverflow.com/a/22470158/19410
[6]: https://medium.freecodecamp.org/how-and-why-you-should-use-python-generators-f6fb56650888
Question: Here is my attempt:

def remove_too_long_reads(bam):

for read in bam:
if read.end - read.start < 60:
yield read

bam = BedTool(input.bam)
bam = BedTool(remove_too_long_reads(bam))

But when I try to use the bam file afterwards like so:

bam.intersect(...

I get an error:

> : Broken pipe The command was:
>
> bedtools intersect -b data/gencode_annotation/tss.bed -a stdin
>
> Things to check: Error in job intersect_data_tss while creating output
> file data/tss/bam/Exp2_9h_PolII.bam. RuleException: KeyError in line
> 31 of
> /local/home/endrebak/code/programmable_epigenetics/rules/tss/intersect_bam_gencode.rules:
> 32 File
> "/local/home/endrebak/code/programmable_epigenetics/rules/tss/intersect_bam_gencode.rules",
> line 31, in __rule_intersect_data_tss File
> "/local/home/endrebak/anaconda3/lib/python3.5/site-packages/pybedtools/bedtool.py",
> line 773, in decorated File
> "/local/home/endrebak/anaconda3/lib/python3.5/site-packages/pybedtools/bedtool.py",
> line 336, in wrapped File
> "/local/home/endrebak/anaconda3/lib/python3.5/site-packages/pybedtools/helpers.py",
> line 373, in call_bedtools

Answer: I have this function in one my scripts - maybe it will help:

```

from pybedtools import BedTool

def filterReadsByLength(inbam, minlength, maxlength):
'''
Takes BedTool which is a bam, and selects intervals that are within the defined lengths.
Input: bam file and min/max lengths
Output: bedTool
'''
print 'Filtering reads by length'
# convert bam to bed
intervals = inbam.bam_to_bed()
filt = intervals.filter(lambda x: len(x) > minlength and len(x) < maxlength).saveas()
return filt

inbam = BedTool(path_to_bam)
filtered = filterReadsByLength(inbam, minlength, maxlength)

```
Question: Dear all,

i am writing to ask please for your suggestions about the appropriate use of CELLRANGER AGGR function in the following context :

we do have 2 experimental replicates for a WT and DISEASE mice :

-- WT_replicate1

-- WT_replicate2

-- DIS_replicate1

-- DIS_replicate2

in order to combine the replicates, shall we just use CELLRANGER AGGR, or shall we just put the fastq files in the same folders i.e.

the fastq files of *WT_replicate1, WT_replicate2* in a folder,

and

the fastq files of *DIS_replicate1, DIS_replicate2* in another folder ?

thanks a lot,

bogdan

Answer: You should not process the FASTQs together, because cell barcodes will clash. So you should process them separately. I would not recommend merging the replicates after processing, information about variance within condition will be useful and important for downstream analysis. If you do aggregate the datasets, you should be aware that the aggr function performs a normalization of the data. In the kallisto|bustools workflow we have a [tutorial][1] for aggregation that allows you to explore other normalizations.

[1]: https://www.kallistobus.tools/aggr_tutorial.html
Question: I have assembled a bacterial genome using spades.
The coverage among the contigs is very heterogeneous, while some have more than 30 fold of coverage, other have less than 5 fold. Also, there are a lot of contigs with a length less than 500 bp.
I'm aware that I must filter contigs having low coverage and small length (they probably are spurious sequencing products), however what would be reasonable criteria for filtering them?
I have found scripts that remove contigs having less than 500 pb and 2 fold coverage. I believe this setting is too permissive, since most of the contigs of my assembly significantly deviate from it.

EDIT:
![Here is the distribution of log10 contig length:][1]

![Here is the distribution of contig coverage:][2]

As we can see, most of contigs have less than 1000 bp, Most of these contigs have less than 5 fold of coverage. My plan is to filter contigs having less than 1000 bp and 5 fold of coverage. After, I am going to check the completeness of the filtered genome sequence using checkm (contrasting to the unfiltered genome).

EDIT 2:
![Length vs Coverage][3]

[1]: https://image.ibb.co/iLswOb/semia4088_log_leng.png
[2]: https://image.ibb.co/i7sMOb/semia4088_cov.png
[3]: https://image.ibb.co/jVC6PG/semia4088_len_vs_coverage.png

Answer: [BlobTools][1] is a great tool to investigate contamination in (expectedly) single-organism assemblies.

[1]: https://github.com/DRL/blobtools