cannin/biostars_qa huggingface dataset:
Question: Hi there,

I have a tumor and a normal BAM file and am preparing to run base recalibration.

I was planning on calling variants on the normal and using that, in addition to dbSNP, as input for recalibration of tumor BAM(s), e.g.:

```
gatk BaseRecalibrator \
-I tumor.bam \
-R hg38.fasta \
--known-sites normal.vcf \
--known-sites dbSNP_hg38.vcf \
-O tumor_recal.table
```
Before producing the normal VCF however, **it's not clear to me whether I should run base recalibration on the normal BAM**. If this is advised, I had planned using dbSNP as the known sites (for normal), e.g.:

```
gatk BaseRecalibrator \
-I normal.bam \
-R hg38.fasta \
--known-sites dbSNP_hg38.vcf \
-O normal_recal.table
```

**Alternatively, I could keep things simple and run base recalibration on both tumor and normal using dbSNP only.**

Is one of these workflows more preferable?
Any clarity here would be much appreciated. Thanks!

Answer: The current "best-practice" is to always do BQSR with the latest ([and largest][1]) dbSNP VCF on all samples - tumor or normal, FFPE or blood, etc. Per discussions in [this post][2] and [this post][3], BQSR can benefit slightly if you provide a "bootstrap of known variants" unique to your samples, either somatic/germline variants found in your tumor/normal. However, you are effectively running your primary analysis pipeline twice (which is overkill), potentially amplifying false-positive variants (from your first-pass variant list), and potentially breaking compatibility of your BAMs with secondary analysis pipelines (e.g. downstream false-positive filters that use BQ).

There are also recent arguments against using BQSR at all, and instead flagging false-positives based on [Base Quality drop-off][4] (at the ends of reads, strand-bias, etc.). You can find an [old Perl script here][5] that implements such BQD filters. I also just found this tweet from Geraldine of GATK saying they're thinking of [dropping BQSR from the best-practices][6] - presumably because the high computational-expense of BQSR is not reasonable when the quality of DNA sequencing has improved so much. I would still recommend BQSR when re-analyzing old FASTQs, or when comparing FASTQs from a mix of different sequencers (e.g. HiSeq and NovaSeq).

[1]: https://www.semanticscholar.org/paper/Effect-of-Database-Size-in-the-Genetic-Variants-Kim-Lee/448b6a85026b3e1a0334d315aebf8ce77269afba
[2]: https://www.biostars.org/p/104852/
[3]: https://sites.google.com/a/broadinstitute.org/legacy-gatk-forum-discussions/2017-06-18-2017-01-18/9713-BQSR-with-MuTect2-use-it-or-not-
[4]: https://support.illumina.com/content/dam/illumina-support/help/Illumina_DRAGEN_Bio_IT_Platform_v3_7_1000000141465/Content/SW/Informatics/Dragen/BaseQualiltyDropoff_fDG.htm
[5]: https://github.com/genome/fpfilter-tool/blob/master/fpfilter.pl#L123-L133
[6]: https://twitter.com/vdageraldine/status/1296181178534440963
Question: Dear all,

I have FASTQ files and on start of my read I have 7 nucleotides tag - I would like to extract reads with this specific tag.

I would like to search in first 15 nucleotides of my reads, if match - extract this read to new fastq files.

Thank you for any ideas or help.

Answer: Probably not the fastest option:

grep -B1 -A2 "^AGATCGG" file.fq | grep -v "^--$" > out.fq

Find lines that begin with "AGATCGG", grab one line before each hit and two lines after each hit, remove lines that are "--".

Note that it's possible (but extremely unlikely) that some quality value line begins like "AGATCGG", in which case the above command would mess up the output file. The likelihood of this is probably very close to zero but if you're processing a file with googolplex lines, maybe it could happen.
Question: I'm struggling with this issue:

./hisat2 -x ../../data/index/hisat_2/iwgsc -U ../../data/mrcv/genes/21dpi_C_R1.fq.gz,../../data/mrcv/genes/21dpi_C_R2.fq.gz,../../data/mrcv/genes/21dpi_C_R3.fq.gz,../../data/mrcv/genes/21dpi_C_R4.fq.gz,../../data/mrcv/genes/21dpi_T_R1.fq.gz,../../data/mrcv/genes/21dpi_T_R2.fq.gz,../../data/mrcv/genes/21dpi_T_R3.fq.gz,../../data/mrcv/genes/21dpi_T_R4.fq.gz -S ../../data/mrcv/genes/alignment/21.sam --un ../../data/mrcv/genes/alignment/21UN.sam -p 6

Error: Read J00113:78:H3JKNBBXX:3:2209:5670:2281 1:N:0:CGATGT has more quality values than read characters.
terminate called after throwing an instance of 'int'
Aborted (core dumped)
(ERR): hisat2-align exited with value 134

Scythe was used to remove adapter contamination and Sickle to remove low-quality reads.

I'm using HISAT2 version 2.1.0 on Ubuntu 18.04. Any ideas what could be wrong?

Answer: I could verify with

gzip -t R2.fq.gz && echo ok || echo bad

that the file was corrupted. For some reason the reads were not in the correct line. Re-downloading the file solved the issue.
Question: Hello Everyone,

I have a bash Script of FastQC and Trimmomatic. I want to know how much resources each of these tool uses. So is there any way or any bash command or tool which helps with this situation.
If possible I would like to change the resources they use to make them perform faster.

Thank you, Nikhil

Answer: You cannot really 'change the resources' they use unless the tool itself has options for the number of threads/memory usage etc.

The only real changes you could make that might marginally speed something up would be to parallelise it with `gnu parallel` (if possible), or maybe change the [`niceness` of the process][1].

As Dariober said, there are not really many options in 'default' tools to monitor resource usage other than at the overall system level once a script is running (e.g. `free`, `htop/top` and `ps` etc). You would have to look at installing some sort of server monitoring solution, but its probably not worth the effort. There are some tools here, but even they are not 'in-script' measures https://www.tecmint.com/command-line-tools-to-monitor-linux-performance/

You can 'profile' your commands by using `time` to figure out where the slow bits are, but unless your code is very inefficient you won't be able to do much to change it.

[1]: https://askubuntu.com/questions/656771/process-niceness-vs-priority
Question: Hi,

I would like to remove human reads from a metagenomics NGS dataset ( fastq ).
Which parameter from **bwa mem** can I change to increase the sensibility ? ( I would like to align more reads )
Penality score ? Gap open ? Clipping penality ?

Answer: Personally, I would use the [RemoveHuman][1] approach from BBMap.

Are you trying to clean Illumina reads? BWA is already pretty sensitive, you will probably increase false positives if you tweak the parameters. If you want increased sensitivity, with (I believe) the same precision, you can try `vg`:

![VG_graph][2]

[1]: http://seqanswers.com/forums/showthread.php?t=42552
[2]: https://cloud.githubusercontent.com/assets/145425/24107859/370225f6-0d8c-11e7-8af4-ab9250971ec8.png
Question: Hello everyone,

I am currently working on a metagenomics analysis using ONT and I have used epi2me wimp for classification of the reads. What I want is to plot the relative abundance of each species present in the sample.

Epi2me outputs a read count assigned to each taxa so I can produce a plot for the proportion of reads per taxa but since each reads has a different length should I normalize by read length ?

Answer: The more relevant normalization would be

1. normalize to one million sequenced reads per sample (if you sequenced more reads, you'll get more hits)
2. normalize to bacterial genome size (if the bacteria has a big genome, you'll get more reads from it compared to one with a small genome)

We cover these in our pipeline Wochenende https://github.com/MHH-RCUG/Wochenende
Question: I have a fasta format file. It looks like. test.fa

>WP1_11_x10
CCATGCGCGGGTTCAATTCCTGTCGTTCGACCC
>WP1_16_x1
ACAAGCGAAGGCTCCTCAACGACGCCTCATCGGATG
>WP1_17_x1
AGAACAAGATTGTTGAAAACTTGAGGA
>WP1_21_x11
CCTGGGATGCGCAAGGAAGCTGAC

I want the below format.The number(10,1,1,11) is in the fasta header .

CCATGCGCGGGTTCAATTCCTGTCGTTCGACCC 10
ACAAGCGAAGGCTCCTCAACGACGCCTCATCGGATG 1
AGAACAAGATTGTTGAAAACTTGAGGA 1
CCTGGGATGCGCAAGGAAGCTGAC 11

Who can help me?Thank you very much!!!

Answer: hai 21.., Please refer to code below written in perl

#! /usr/bin/perl -w

open(F1, '<', "$ARGV[0]") or die;

while (<F1>){

if ($.%2 != 0) {

@val = split(/\_/);
$cnt = $val[1];
}
else{ chomp($_); print "$_ $cnt \n"; }
}

This code was not written in a defensive way, this should work normally with out any problems.
Question: I have some bam files that I want to split. The reference genome has 22 chromosomes and a bunch of unplaced scaffolds. I can easily split the bam files chromosome by chromososme

samtools view my.bam chromsosmeI -b -o my_chromosomeI.bam

I can do the same *scaffold by scaffold*, but I would like to do extract all scaffolds with some sort of wild card (there are a *lot* of them) and have them all in the same bam output. The wild card part is the one I'm having issues.

Answer: samtools idxstats in.bam | cut -f1 | grep | xargs samtools view -o scaffolds.bam in.bam
Question: Hi all,

I'm trying to build a custom reference panel for use with IMPUTE2. I have a phased .vcf from all my participants, and I was using [the first script here][1], modified slightly to work with my vcf/HG38, to build the reference panel. The exact command I run to generate the files for chromosome 1 (as an example) is:

./vcf2impute_legend_haps.pl -vcf my_vcf.vcf.gz \
-leghap output_prefix \
-chr chr1 \
-hg38

The -hg38 thing is just an argument I added to the script, and can be ignored.

This generates .hap.gz files, .legend.gz files, and .sample_list files, one for each chromosome. What I'm missing is a genetic map file. I don't what these are, and I can't find anything on it in the documentation. Looking at the 1000 Genomes HG37 reference panel files provided [here][2], it's not clear to me what these mapping files are, so I can't replicate the steps to create them.

Any advice would be appreciated.

[1]: https://mathgen.stats.ox.ac.uk/impute/impute_v2.html#scripts
[2]: https://mathgen.stats.ox.ac.uk/impute/impute_v2.html#reference

Answer: According to the documentation you linked, genetic maps are recombination maps, measured in centiMorgans.

If you can lift-over the genomic position of the markers from the 1000 Genomes HG37 reference panel, you can use the same genetic map position (the **Genetic_Map(cM)** column), and recalculate the recombination rate between current position and next position in map (the **COMBINED_rate(cM/Mb)** column).

However, are you sure there isn't an updated version of these files?
Question: I have a question regarding stranded RNA-Seq from Illumina TruSeq. I'm comparing some stranded RNA-Seq runs from SOLID and Illumina. I have found that while SOLID stranded protocol is perfect (literally no reads aligning to antisense strand), I found quite a lot reads in antisense strand in Illumina sample. Still sense strand reads are enriched (~20 times more reads in sense strand).

The reads from antisense strand are aligned uniquely (MAPQ>100, mapped with STAR), but I don't think it's real expression from antisense strand, as the exon-intron structure of antisense reads is similar to exon-intro structure from sense. But I see huge enrichment 10-20x) in Illumina

Is TruSeq designed to enrich for sense strand or is it suppose to be perfectly stranded protocol as SOLID was?

Is it possible that something went wrong during library prep?

Answer: I have found highly variable sense-strand read enrichment (6.39-27.05; median: 9.85) in our RNA-Seq stranded data (6 time points with 3 RNA fractions in colours). I have ignored exons having any overlap with known antisense exons.

Overall, red RNA fraction tends to have higher sense-strand read enrichment that blue/orange. **Any clue why there is such high variability in enrichment between time points / RNA fractions?**

I have obtained slightly lower enrichment for genes (5.98-24.07; median: 8.74). As these data is from zebrafish, for which the annotation is far from being perfect, I suspect many antisense genes / exons are simply not annotated.

More about methodology [here][1].

[1]: http://bioinfoexpert.com/?p=447
Question: Hi,

I want to convert my files (xx.bed, xx.bim, xx.fam) to vcf format. I am following the code given at: http://pngu.mgh.harvard.edu/~purcell/plink2/data.html#recode

Here is my command line output:

```
============================
node00:[exomes] % ls
exomes.bed exomes.bim exomes.fam

node00:[exomes] % plink --bfile exomes --recode vcf-iid --out ./new_vcf

@----------------------------------------------------------@
| PLINK! | v1.07 | 10/Aug/2009 |
|----------------------------------------------------------|
| (C) 2009 Shaun Purcell, GNU General Public License, v2 |
|----------------------------------------------------------|
| For documentation, citation & bug-report instructions: |
| http://pngu.mgh.harvard.edu/purcell/plink/ |
@----------------------------------------------------------@

Web-based version check ( --noweb to skip )
Recent cached web-check found...Problem connecting to web

Writing this text to log file [ ./new_vcf.log ]
Analysis started: Tue Apr 21 15:02:41 2015

** Unused command line option: vcf-iid

ERROR: Problem parsing the command line arguments.
==========================================
```

I also tried giving the options 'vcf' and 'vcf-fid' but got the same results.

What am I doing wrong? Or is there any other tool that I can use to achieve this?

Thanks

Answer: You are using PLINK 1.07, and this command requires PLINK 1.9.
Question: Hello Biostars,
Can anyone tell me how to prepare input data set for GSEA after Differential Gene Expression Analysis by DESeq2?
How will I rank the genes? Should I rank based on log2FC or Adjusted P value? Is there any way to generate a GSEA ready data directly from DESeq2?.
I was using topGo for gene ontology enrichment analysis before and recently came across GSEA. Which one is better GO enrichment analysis or GSEA? Even after going through the papers I couldn't find a significant difference between above two.

Thank you

Answer: Here's a link to an answer I wrote a few years ago for using the gene set testing package *goseq* following *DESeq2*:

https://support.bioconductor.org/p/64811/#64815

I'm not sure what kind of input GSEA takes. I also like the methods behind *ROAST* and *CAMERA* from the limma package, but I haven't yet worked on integrating with those methods. For those two, you would need to run a limma analysis upstream.
Question: Hi,

I have imputed genotypes using the IMPUTE2 software. Now I need to convert my IMPUTE2 formatted genotypes into VCF. I have used the qctool to do so. However, the resulting VCF file has lost the phasing information. Is there any tool that allows me to convert IMPUTE2 to VCF while keeping the phasing information?

Thanks

Answer: I ended up using shapeit (version v2.837) -convert function to convert IMPUTE2 to VCF while keeping phase information.
Question: Hello:

I explore between studies about associated histone mark in enhancer locations and the result is very diverse and confusing. I am looking for enhancer STAT1. I would like to know the histone marks associated with it?

thanks.

Answer: H3K4me and H3k27ac are histone marks for enhancers. p300 (a histone acetyl transferase) is also found at enhancers.
Question: I have a fasta file with hundreds of sequences and their respective headers. The headers (all of them) are in the format

>ABCD [id_123] (gene_XYZ) [protein_ijk] [protein_id=qqq] [123..899]
.......seqeunce............

>EFGH [id_999] (gene_PQR) [protein_tre] [protein_id=trs] [573..789]
......seqeunce............

and so on.....

For the header every info in parenthesis are continuous and are only separated by a single space each (just as written above).
All I want to do is retain "ABCD" (the very first info) in the header corresponding to every sequence . I want to loop through all the headers that are present in the file and return something like this :

>ABCD
.....sequence.....
>EFGH
.....sequence.......

and so on......Any help is most appreciated and i am working with BASH and perl.

Thank and regards!

Answer: In bash: `cut -d ' ' -f1 your_file.fa > new_file.fa`
Question: I have a list of protein GI and would like to get accession number of the genome (DBSOURCE in genbank file) using Eutilities.

What would be the easiest way to get genome accession numbers for a list of protein GI?

Answer: An alternative to perl e-util is to use the Unix e-utils. Following one liner does the job!
Note that following command would work for accession number as well as GIs as -id parameter in elink command.

elink -db protein -id 817524604 -target nuccore|efetch -format acc
Question: Hello,

I have this bam file that originally contained paired sequences, after converting it into fastq sequences I got 2 different files. My idea is to make an alignment with a chromosome sequence, and from there obtain a VCF. My question is if the sequence that I will use will matter in my future results?, there is a way to know which sequence should I use?

I got the sequences using samtools:

bam2fastq -o blah_aligned_R#.fastq --no-unaligned Bamfilename.bam

EDIT:

Thank you very much for the comments, I will explain myself better (and sorry for any misunderstanding, but this is the first time that I am working in bioinformatics).

I am wanting to see the occurrence of KRAS-variant mutation (as described by Ratner et al, in http://www.nature.com/onc/journal/v31/n42/pdf/onc2011539a.pdf) in chromosome 12 in ovarian cancer cell lines. For this, I have downloaded some bam files from the nci60-project (http://watson.nci.nih.gov/projects/nci60/wes/BAMS/). These bam files are presented as paired sequences, so after extracting chromosome 12 information from them (using the bai files that are in the same webpage as index) by the code:

$ samtools view IGR-OV1_reord_mdups_ralgn_fmate_recal.bam chr12>chr12.bam

Then, I sorted the bam file, and instead of using the original code that I posted at the beginning, I converted them to fastq using:

```
$ bedtools bamtofastq -i chr12.qsort.bam \
-fq chr12.end1.fq \
-fq2 chr12.end2.fq
```

So now, I have two sequences in fastq format per each original bam file.

My following steps, that I have in mind, is to align these fastq files with the human assembly 37, convert then again to bam files, and make the variant calling with mpileup so from there start with the annotation using annovar or something alike. So my questions in fact are:

1. Is my sort of pipeline acceptable?
2. Should I use both sequences for the variant calling?

Thank you in advance for your help and comments

Answer: The steps sound fine, and as others have suggested, using both read1 and read2 for alignment is more efficient. The paired mapping allows for better read placement, especially over repetitive or homologous regions.

Am not sure why you dont directly use the aligned BAM file to call variants on chr12 and then annotate them or intersect for any falling in KRAS?

If it is just a single position in KRAS, seems the paper you ref its a 3' UTR mutation, you can simply genotype that coordinate (chromosome-position) on all the BAM files and get a quick answer if the mutation is present/absent.

Hope this helps.
Question: Aloha!

I have a fasta file that looks like the following:

>FFSA34B_100_M7_ID10014
ATCTAACAATGTTGCTCATGCAGGCCCTGCAGTAGATTTAACCATTCTATCCCTTCACCTAGCAGGTGTATCCTCCTTAATAGGAGCCATCAATTTTACAACTACTATTGCTAACAGACGTTTAGAAGGTATACCTACAGAAAAAATACCCTTATTTATT
>ART01B_100_M7_ID10021
TGTCTTCTAATATTTCCCATGGAGGGGCCTCTGTTGATATGGCTACTTTTTCTCTCCACCTGGCTGGTGTGTCCTCTATTTTGGGTGCTATCAATTTTATAACTACTATTATTAATATGCGGAGTAATTCCTTGTCTTTTGTGCGGGTTCCTTTGTTTGT
>PHR01B_100_M7_ID100236
TAAGAGGAAATGTAGCCCATAGTGGACCTGCAGTAGATTTAGGAATTTTTTCTTTACATTTAGCTGGGGTTTCTTCTATTTTAGGGGCTATTAACTTTATTACTACAATTATTAGTATACGGGAGAATATGGGGATAGAAAAAGTTACTTTATTTAGTTG

I have a text file that looks like the following:

FFSA34B_100_M7_ID10014 2
ART01B_100_M7_ID10021 454
PHR01B_100_M7_ID100236 43

And I would like to append the second column of the text file to the matching fasta header to produce the following output:

>FFSA34B_100_M7_ID10014_2
ATCTAACAATGTTGCTCATGCAGGCCCTGCAGTAGATTTAACCATTCTATCCCTTCACCTAGCAGGTGTATCCTCCTTAATAGGAGCCATCAATTTTACAACTACTATTGCTAACAGACGTTTAGAAGGTATACCTACAGAAAAAATACCCTTATTTATT
>ART01B_100_M7_ID10021_454
TGTCTTCTAATATTTCCCATGGAGGGGCCTCTGTTGATATGGCTACTTTTTCTCTCCACCTGGCTGGTGTGTCCTCTATTTTGGGTGCTATCAATTTTATAACTACTATTATTAATATGCGGAGTAATTCCTTGTCTTTTGTGCGGGTTCCTTTGTTTGT
>PHR01B_100_M7_ID100236_43
TAAGAGGAAATGTAGCCCATAGTGGACCTGCAGTAGATTTAGGAATTTTTTCTTTACATTTAGCTGGGGTTTCTTCTATTTTAGGGGCTATTAACTTTATTACTACAATTATTAGTATACGGGAGAATATGGGGATAGAAAAAGTTACTTTATTTAGTTG

The fasta headers and the text file in the real example will not be in the same order so the command will have to match the fasta header to the first column of the text file and then append the second column value to the fasta header. Anyone know of a why to easily do this?

Thanks in advance!

Answer: [linearize][1], sort and join (assuming delimiter is tab)

$ join -t $'\t' -1 1 -2 1 \
<(awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' < input.fa | cut -c 2- | sort -t $'\t' -k1,1) \
<(sort -t $'\t' -k1,1 list.txt) |\
awk -F '\t' '{printf(">%s_%s\n%s\n",$1,$3,$2);}'

>ART01B_100_M7_ID10021_454
TGTCTTCTAATATTTCCCATGGAGGGGCCTCTGTTGATATGGCTACTTTTTCTCTCCACCTGGCTGGTGTGTCCTCTATTTTGGGTGCTATCAATTTTATAACTACTATTATTAATATGCGGAGTAATTCCTTGTCTTTTGTGCGGGTTCCTTTGTTTGT
>FFSA34B_100_M7_ID10014_2
ATCTAACAATGTTGCTCATGCAGGCCCTGCAGTAGATTTAACCATTCTATCCCTTCACCTAGCAGGTGTATCCTCCTTAATAGGAGCCATCAATTTTACAACTACTATTGCTAACAGACGTTTAGAAGGTATACCTACAGAAAAAATACCCTTATTTATT
>PHR01B_100_M7_ID100236_43
TAAGAGGAAATGTAGCCCATAGTGGACCTGCAGTAGATTTAGGAATTTTTTCTTTACATTTAGCTGGGGTTTCTTCTATTTTAGGGGCTATTAACTTTATTACTACAATTATTAGTATACGGGAGAATATGGGGATAGAAAAAGTTACTTTATTTAGTTG

[1]: https://gist.github.com/lindenb/2c0d4e11fd8a96d4c345
Question: Is it possible to skip the estimateSizeFactor step in Deseq2 and define your own size factors? If I just replace the size factor vector with my own values, will the down-stream DE steps use these values? I am pretty new with R so it would be very helpful if I can see some example code on how to do this.

Answer: Hi,

Short answer : yes and yes !

As for the example code, see my answer [here][1].

[1]: https://www.biostars.org/p/159472/#159511

or more to the point :

sizeFactors(cds) = c(my_Values)
Question: Hello,

I have been using VarScan and MuTect to call somatic variants. Lately, of curiosity, I also tried to call variants separately for the tumor sample (vs hg19) and normal sample (also vs hg19) using GATK, and compared the outputs, which is - I thought - theoretically also the somatic variants. However, the number of variants found this way is way more than the one found by VarScan or MuTect.

My questions is: What's wrong with doing so? Is that because some systematic errors in the algorithm somehow got doubled when calling variants separately?

Thank you very much.

Answer: There are posts that hint at a lot of this if you do a search on here. The short answer is that somatic mutation calling is a different beast. To give just one example, single sample (germline) calling, you can assume that the frequency of SNPs will generally be 50% or 100%. In tumors, you have to take purity, ploidy, and tumor heterogeneity into account and SNPs may occur across a wide range of VAFs.

We get into all of this a little bit in our recent paper here: [Optimizing Cancer Genome Sequencing and Analysis](http://www.cell.com/cell-systems/abstract/S2405-4712(15)00113-1)
Question: Hello
I am trying to build my index to HISAT2 using the following command

hisat2-build rn6_ercc/genome.1.ht2,genome.2.ht2,genome.3.ht2,genome.4.ht2,genome.5.ht2,genome.6.ht2,genome.7.ht2,genome.8.ht2 ht2_base

I got the following errors

>Error: could not open genome.2.ht2

and

>Error: Encountered internal HISAT2 exception (#1)

Kindly help

Answer: You should build the index with something like:

hisat2-build genome.fa index_name

Where genome.fa is the fasta file of the genome you feed hisat2-build, and index_name will be the base name of the index files (*.ht2).