cannin/biostars_qa huggingface dataset:
Question: I have a SAM file with alignments and for each entry alignment, I want to reconstruct the alignment between the reference and the read based on the CIGAR and MD strings. It seems like this should be possible, but this one example bothers me:

SRR037452.3355 0 ENSG00000266658|ENST00000607521 3523 255 16M1I18M * 0 0 CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC IIIIIIIIIIIIIII:III/IIII=+IGC,I"/I. NH:i:1 HI:i:1 NM:i:2 MD:Z:33G0

Here, the CIGAR string has an insertion to the reference which messes up the MD string indexing. According to MD string, the read should be only 34 bases long (r=35). My guess is that the alignment is actually this: 16 matches, 1 base that is present in the read, not present in the reference, 17 matches, one mismatch where the read has a "G" and the reference has something else. Is that correct? Are there CIGAR/MD string combinations where reconstructing the alignment would be impossible (i.e. either of the strings is ambiguous)?

Answer: Your guess is correct, the alignment has 16 matches, 1 unspecified insertion, followed by 17 matches and one mismatch ('G' on the reference which is substituted with a 'C' in the read sequence).

If the MD tag is present, you should always be able to reconstruct the alignment from the MD tag and the CIGAR, as long as you don't care about the sequence of any mismatched or inserted bases (which you can only get from the read sequence).
Question: Hello, I have phased plink files of a population. I would like to obtain all SNPs that are in high linkage disequilibrium (LD) for a set of 10 SNPs. I know how to do it according to a r2 threshold, but I need to set this threshold according to D'

For instance, this script will output all variants in high LD (r2 > 0.6) for a list of 10 SNPs within a range of 1000 Kb.

plink \
--bfile myfile \
--r2 --ld-snp-list list_10_snps.txt \
--ld-window-kb 1000 --ld-window-r2 0.6

Thanks very much

Answer: This isn't directly supported by plink, but you can use "--r2 dprime" to include D' in the main LD report, and then follow up with e.g. an awk one-liner to filter on D'.
Question: How to embed known cell information into Seurat object? Using AddMetaData?

Is this possible?

I read this documentation: https://www.rdocumentation.org/packages/Seurat/versions/3.1.4/topics/AddMetaData

However I need some clarification on how to go about doing this, please?

I have cell labels that are all aligned with the expression matrix such as this:

"12wks_fetal pancreas cell_acinar"	"19wks_fetal pancreas cell_beta"	"12wks_fetal pancreas cell_ductal"	"12wks_fetal pancreas cell_beta"	"12wks_fetal pancreas cell_acinar"	"22wks_fetal pancreas cell_acinar"	"12wks_fetal pancreas cell_acinar"	"12wks_fetal pancreas cell_alpha"	"12wks_fetal pancreas cell_ductal"	"12wks_fetal pancreas cell_acinar"	"14wks_fetal pancreas cell_endocrine.progenitor..."

I used these cell labels as cell names (column names) for the expression matrix. However that's not useful for seeing which cell is from what time point 12wks, 19wks, 14 wks, etc. on the UMAP for example.

Would there be a way to label the cells using the AddMetaData function so that when I visualize the data using UMAP, they are color coded based on the time of collection (ie. 12wks, 19wks, 22wks, etc.) or would there be a way I could also filter by both time of collection and/or cell type?

Is this possible in Seurat? How would I approach this?

I would really appreciate anyone's help.

Very Respectfully,
Pratik

Answer: You didn't think to try AddMetaData?

Yes, if you have one column with barcode names an another with any kind of information, you can add a new column of metadata. You can also make a new column of metadata by concatenating two columns you already have. You can use group.by in DimPlot to color the cells by any column in the metadata, and you can subset by metadata columns after setting that column with Ident.
Question: Probably a basic question but it's difficult to find information on the subject...

When viewing a GTF and BED file in IGV, there seems to be differences in colour (BED track can't be recoloured, just shows as black, GTF can be recoloured). Does anyone know the reason for this? I thought that the GTF and BED formats contained predominantly the same information just in different formats.

Thanks!

Answer: A GTF file is more akin to a large group of BED files grouped according to some feature or features of relevance, often with multiple levels of hierarchical metadata. In the classical case of genes, GTF specifies a standard way of relating exons to transcripts and genes. There's no standard way of doing that with BED files. BED files are useful for cases where the file itself contains only a single type of information (e.g. ChIP-seq peaks), where you could instead combine multiple of these files in a GFF (e.g., annotating peaks from different histones pulled down with multiple antibodies each). You can do some of this by using the `name` field in the BED format, but there's still no good way to represent hierarchical.
Question: hi,

I have a list of miRNAs IDs from different organisms how I can take sequence of mature miRNAs for all my IDs?

for example

lgi-mir-34

mdo-mir-193

mdo-mir-21

mdo-mir-302d

mml-mir-18

Answer: another way was mirdeep2

fastaselect.pl reads.fa reads_select.ids > reads_select.fa

This script only prints out the fasta entries that match an id in the id file.
Question: Hello friends,
I recently analysed a RNA-Seq data and yielded DEGs using DESeq2.
Everything was fine until I figured out that on applying P value 0.05, i am left with 1,700 DEGs whereas on considering and applying Padj 0.05, I am left with only 127 DEGs??? isn't this surprising ? How can there be such a huge difference in number of significant data on just considering adj P value?

Thank you to anyone who can afford any little energy and time to guide me. I am stuck!!

Answer: Its completely expected. Under classic multiple testing correction you'd expect your adjusted P-values to be 20,000 times larger than you non-adjusted values. Obviously its not as bad as this because we don't use Bonferroni in RNAseq, but the multiple testing burden is still high.

Another way of looking at it is that under a standard hypothesis test, we would expect 5% of tests to give a false positive, thats around 1000 genes, so we would expect the number of genes passing 5% threshold to be at least 1000 less for an adjusted p-value.
Question: Hi all,

I am very new in this filed and really appriciate whoever can help me with this question. I was wondering if you how I can filter my bfiles( BIM, BED, FAM) for a subset of my SNPs. I have my BIM , BED, FAM files. But I do not need all SNPs available so I filterd my BIM file by AWK commands but I do not know how to match my BED file with the new BIM file created. Please tell me how should I do it? thank you very much in advance.

Answer: you should create a list of your snps and use plink's filtering option. I believe it goes something like :

plink --bfile --remove [snps.txt] --out

this should also subset your BED file
Question: **hello,
when i run gatk, the error always occur, like this "java.lang.IllegalArgumentException: samples cannot be empty"
is there mistake in my input file?
thank you for your help !**

**my bam file as following:**
HWI-EAS418:3:37:1070:1462	83	chr20	46689301	255	50M	=	46687222	-2129	CAGCTCCAGGCCGCTCAAGAAGCGGCTGCTCCGCTCCCGGGCTGCGGCCA 0:@2%,6.=:.4=,4+6;B79>=BB=2>7=:9=4(599@BB9BB>BB??A	NH:i:1	HI:i:1	AS:i:99	nM:i:0

**the script is following:**
java -jar /home/H/mutation/gatk/gatk-package-4.1.4.1-local.jar HaplotypeCaller -R /home/H/mutation/ref/ctat_genome_lib_build_dir/ref_genome.fa \
-I /home/H/mutation/ctat-mutations-master/testing/__misc_data/Aligned.sortedByCoord.out.GRCh38.bam \
--recover-dangling-heads true \
--dont-use-soft-clipped-bases \
-stand-call-conf 20.0 -O test.vcf

**the logs are followinig:**
Runtime.totalMemory()=2260729856
java.lang.IllegalArgumentException: samples cannot be empty
at org.broadinstitute.hellbender.utils.Utils.validateArg(Utils.java:725)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.ReferenceConfidenceModel.(ReferenceConfidenceModel.java:116)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerEngine.initialize(HaplotypeCallerEngine.java:221)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerEngine.(HaplotypeCallerEngine.java:164)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller.onTraversalStart(HaplotypeCaller.java:224)
at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1046)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:139)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:191)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:210)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:163)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:206)
at org.broadinstitute.hellbender.Main.main(Main.java:292)

Answer: I had the same problem, I had to add the read group @RG tag to the .bam file using

`samtools addreplacerg -r '@RG\tID:samplename\tSM:samplename' input.bam -o output.bam`
Question: I m looking for a source code that can take multiple sequence alignment and return back the Sum-of-pairs or weight sum-of-pairs score based on any scoring scheme.. any one can help..

Answer: when I moved on from CLUSTAL, I found MAFFT to be a great (and amazingly faster) alternative for multiple sequence alignment that even allowed to play around with scoring matrices. I know it takes sum-of-pairs into account when setting up a run, but what I don't know is if it provides them as output. the source is available and published, so I think it would be worth it giving it a try.
Question: Hi all,

I have recently tried to estimate runs of homozygosity (ROH) from my vcf file by using plink 1.9.

I ran following code to generate binary files that plink required:

plink --vcf myfile.vcf --make-bed --out out_name --no-sex --no-parents --no-fid --no-pheno --allow-extra-chr

This vcf file only contains one individual and about 3 million SNPs.

I used --allow-extra-chr here because I mapped my sequences to a drift genome.

Then, I used following code (with default parameters), trying to estimate ROH of my sample:

plink -bfile out_name --homozyg --allow-extra-chr

The result gave me 0 ROH, and only header produced in .hom file.

I also tried different parameters with different SNP windows and criterions, such like:

plink -bfile out_name --homozyg --homozyg-window-snp 50 --homozyg-snp 50 --homozyg-window-missing 3 --homozyg-kb 100 --homozyg-density 1000 --allow-extra-chr

However, all the results were the same, that showed :

PLINK v1.90b6.12 64-bit (28 Oct 2019) www.cog-genomics.org/plink/1.9/
(C) 2005-2019 Shaun Purcell, Christopher Chang GNU General Public License v3
Logging to plink.log.
Options in effect:

--allow-extra-chr
--bfile out_name
--homozyg
--homozyg-density 1000
--homozyg-kb 100
--homozyg-snp 50
--homozyg-window-missing 3
--homozyg-window-snp 50

515905 MB RAM detected; reserving 257952 MB for main workspace.
3708761 variants loaded from .bim file.
1 person (0 males, 0 females, 1 ambiguous) loaded from .fam.
Ambiguous sex ID written to plink.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1 founder and 0 nonfounders present.
Calculating allele frequencies... done.
3708761 variants and 1 person pass filters and QC.
Note: No phenotypes present.
--homozyg: Scan complete, found 0 ROH.

Results saved to plink.hom + plink.hom.indiv + plink.hom.summary .

Does anyone have idea why this happened to my files?

Will be grateful for any suggestions.

Answer: This is unsurprising if your VCF contains only positions where your sample *differs* from the reference genome; practically all the ROHs will span regions which are excluded by your VCF. You may need to reconstruct your VCF in a way that includes homozygous-REF calls in an unbiased manner.
Question: Hi, guys please i need some help because i am not confident with linux
first i have managed to installed ubuntu on a vitual box on my windows 7 and now
i want to install vcf tools and make exeutable files
please how can i do this.
thanks

Answer: This is from the installation instruction in the read me:

> Compiling:
> To compile and install VCFtools, you should type 'make' in the vcftools folder. The perl scripts and cpp executable will be installed in the /vcftools_(version_num)/bin/ folder. It is recommended that you add this folder to you PATH.

So, go to the folder where you have extracted vcftools by typing:

cd [path-to-vcf-tools]
make

Then you can add it to your path, if you like, (so that you can run the tools by typing there name on the command line, rather then referencing them exactly) by following the instructions [here](http://askubuntu.com/questions/60218/how-to-add-a-directory-to-my-path).

Alternative way:
If you have the Ubuntu universe repositories enabled you can install it from there (of course assuming that you are running Ubuntu). Just open a terminal and type:

sudo apt-get install vcftools

Hope you get it working.

Edit 1: Fixed formatting.

Edit 2: Added easier way to install vcftools in from the Ubuntu universe repositories
Question: I am trying to extract maximum features for my machine learning input, and I need PV4, BQ, G3, NS to be added to my VCF files for all the SNP calling
I used the sametools and bcftools as below to extract some of the information:

samtools mpileup --skip-indels -m 1 -E --BCF --output-tags DP,AD,ADF,ADR,SP,INFO/AD,INFO/ADF,INFO/ADR -f test/reference/Aradu.fa test/bam/A72.bam | bcftools call -mv --skip-variants indels --multiallelic-caller --variants-only |bcftools +fill-tags> test/test/At.vcf

The result of my code has

> Aradu.A01	1345	.	T	G	5.57134	.	DP=6;ADF=0,5;ADR=1,0;AD=1,5;VDB=0.0340507;SGB=-0.590765;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;AC=2;AN=2;DP4=0,1,5,0;MQ=9;NS=1;AF=1;MAF=0;AC_Het=0;AC_Hom=2;AC_Hemi=0;HWE=1	GT:PL:DP:SP:ADF:ADR:AD	1/1:31,10,0:6:0:0,5:1,0:1,5

How can I include any of these PV4, BQ, G3, NS to the file

Thanks

Answer: Hi Sahar,

I get **PV4** (*P-values for strand bias, baseQ bias, mapQ bias and tail distance bias*) and **G3** (*ML estimate of genotype frequencies*) automatically when I align my data with `bwa mem` and then call variants with `samtools mpileup` piped into `bcftools call` (latest versions).

**NS** (*Number of Samples With Data*) may be a tag that was used a lot in the past but that has been more or less replaced. You can obtains similar information by looking at the **AC** (*allele count in genotypes, for each ALT allele, in the same order as listed*), **AF** (*allele frequency for each ALT allele in the same order as listed: use this when estimated from primary data, not called genotypes*), and **AN** (*total number of alleles in called genotypes*) tags.

The same may be true for **BQ** (*RMS base quality at this position*). However, using `samtools mpileup`, you can ensure a minimal base quality on variant bases with the `--min-BQ` comman-lime parameter.
Question: Good afternoon,
I'm working with gatk and more specifically with the Haplotypecaller function in order to create the GVCF files and I've seen that it takes quite a lot of time. I'm a bit in a rush and I would like to speed up the process, but I did not find useful information about it. Here is my code:

java -jar ~/softwares/GATKK/gatk/gatk-package-4.1.7.0-local.jar HaplotypeCaller --reference Pmuralis_1.0.fa --input run2_mergeandaligned.bam --output run2_4096_mergeandaligned.g.vcf -ERC GVCF

The only thing that seems to improve a bit the process is if I add `-Xmx4096m` to the beginning, like this:

java -Xmx4096m -jar ~/softwares/GATKK/gatk/gatk-package-4.1.7.0-local.jar HaplotypeCaller --reference Pmuralis_1.0.fa --input run2_mergeandaligned.bam --output run2_4096_mergeandaligned.g.vcf -ERC GVCF

Another thing I noticed is this message:

20:58:54.759 INFO IntelPairHmm - Available threads: 20
20:58:54.759 INFO IntelPairHmm - Requested threads: 4

As in my server I have 20 cores but the process is just taking 4. I think I solved it adding `-native-pair-hmm-threads 20` but it didn't speed up the process...
Let's see if somebody knows about java and can help me!

Thank you very much!

Answer: > . I'm a bit in a rush and I would like to speed up the process, but I did not find useful information about it.

split by chromosome using option `-L` and run in parallel.
Question: Hi,

I am fairly new to bioinformatics. I have the genome sequences of two closely related bacteria (or, more accurately, both genome sequences are publicly available). One bacterium grows much faster under a particular set of growth conditions than the other. I'm trying to determine what genes may be responsible for this improved growth. How do I compare these two genomes to find out what genes are different between the two? And is there a way to find out what the functions (e.g. GO terms) for the different genes are or will I just get a list of gene names as output?

Thanks for any help you can offer.

Answer: Fast option is to look up DB that store informations about microbial genomes, i.e. [IMG][2]

A more detailed analysis would require a bit more time and some command line options.
If your strains are related (i.e. same Genus or Species), I would say that the first thing to do is to calculate the genes that these two bacteria shares and the one that are unique. This means to identify orthologous genes. There are plenty of programs to do that out there. You can take a look at this tools: https://omictools.com/pangenomics-category. They are intended to study pan-genomes, but most of them will perfectly suit your use. I personally use GET_HOMOLOGUES, but I read that [OrthoFINDER][1] is really good and reliable.

After this analysis, you will have a list of genes that are shared between the two bacteria and a list of genes that are unique. Also, you will see which genes are present in duplicate (paralogues) in your genomes.

However, I am not sure that this analysis will answer your question. Especially if the strain are extremely closely related, you might need to do a more detailed analysis considering gene variants, i.e. SNPs. You might then see that some shared genes show differences in the two strains. You might find that these differences will influence the protein sequences as well (they are non synonymous mutations), indicating that the two strains have different version of shared proteins.

Hope this helps.

[1]: http://www.stevekellylab.com/software/orthofinder
[2]: https://img.jgi.doe.gov/
Question: Hi,

I'm running bam-readcount (commit `6c3f3ae901`) on a few hundred bam files against a single reference fasta file (designated with `-f`). However, when I look at the output for any of the files, all of the reference bases, for any position, are N.

NODE_14_length_46_cov_1.239130 52 N 1 =:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 A:1:2.00:36.00:2.00:1:0:0.00:0.02:0.00:0:0.00:0.00:0.00 C:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 G:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 T:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00 N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00

Is this a bug? The base in the reference sequence at this position is A.

>NODE_14_length_46_cov_1.239130
AGCTAACTGAGTTTATCACACTCAGTTAATGTCCATTTCACTTCACACATAACCTTACAG
ATCGGAAGATCTCGTA

Thanks!

Answer: I see this if I run on the entire BAM file without specifying a list of positions or a region. This is a bug. I pushed a very quick fix out that should solve this issue. It's 4b6479a42d002d855eda6a45bca097756d493cdb. Does this fix the issue?
Question: I'm new to the Bioinformatics field and am still learning how to use all the tools and such, but I've been stumped on something that I feel should be very easy for about a week.

Basically, I want to compare the reads in a BAM file to a FASTA reference sequence and get the changes. I can easily view them with **samtools tview** but I need to have it in an Excel file for manipulation/concatenation. I can't seem to figure out why **mpileup** won't give me a straight list of what reads are mutated, unless that's the wrong tool for the job. I'm pretty good at programming, so if the solution is something I need to code for myself I could do it but since the formats are different than what I'm used to, it's been a major headache.

Thanks

Answer: you need to use bcftools to generate a tabular VCF file:

* http://samtools.sourceforge.net/mpileup.shtml (old version)
* http://www.htslib.org/workflow/#mapping_to_variant (new version)

PS: no, you don't want to have it in an Excel file

PS2: no, you really don't want to have it in an Excel file

PS3: Every time you want to have it in an Excel file, God kills a kitten

PS4:
Question: Hi,

I have RNA-Seq data from patients, and I want to find novel transcripts that are differentially expressed in the treatment condition. For this, I followed the methodology given in this Galaxy tutorial - https://galaxyproject.org/tutorials/nt_rnaseq/. Here they use HISAT for mapping, followed by Stringtie for transcript reconstruction (using Stringtie allows them to find novel transcripts), and then use FeatureCounts for counting the number of reads per transcript. Then they use DESeq2 for differential expression analysis.

My question is that since Stringtie itself also gives the count (in terms of coverage, FPKM and TPM) for each transcript that it constructs, then why should I use FeatureCounts for these constructed transcripts?

Answer: There is absolutely no reason to do that. Actually you should rather use StringTie's quantification since that is more accurate than featureCounts - you can read more about such considerations [here][1]. Acutally DEXSeq directly supports analysis of StringTie data via tximport::tximport() as described [here][2] - although they actually don't mention StringTie in the DESeq2 vignette tximport does support StrigTie as describe [here][3].

[1]: https://bioconductor.org/packages/devel/bioc/vignettes/IsoformSwitchAnalyzeR/inst/doc/IsoformSwitchAnalyzeR.html#what-quantification-tools-should-i-use
[2]: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#transcript-abundance-files-and-tximport-input
[3]: http://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html#stringtie
Question: I have a freebayes genotyped multisample vcf file. I want to sort the names of the samples in alphabetical order to make my life a bit easy with post variant calling analysis. I have tried Picared SortVcf which work fine which works fine on a small file but failed on a large file (~4 Gb). If you suggest ways to sort a large multi-sample file without distorting the variants information, it would a great help.

Answer: Hey, try this:

$ bcftools query -l input.vcf | sort > samples.txt
$ bcftools view -S samples.txt input.vcf > output.vcf

If not already done, I would also suggest to use `bcf` instead of `vcf` or `vcf.gz`. This really improves speed when working with `bcftools` on large datasets.

fin swimmer