Hi all,
I am a 1st year MS student and I am required to compare the data I have with the healthy controls, then find p value and pot the Manhattan plot to find out which mutation is abundant for the condition.
The data I have been given has rsid and gene names, EXAC all and EXAC 1000
I want to know:
how do I compare the entire data I have with the data of the healthy control?
How do I then find the p value?
How do I then plot Manhattan Plot?
Please guide me.
Dear jbagaria1206,
The typical program for performing statistical analyses on genetic data is PLINK. To run PLINK, your starting data should generally be VCF, BCF, or some other standard format. Can you clarify the format of your datasets?
It is possible to run statistical analyses in R, too. The most basic association test is just a Chi-squared test, after all, which I show here in R code: A: SNP dataset and Z Score
I also have a Bioconductor R package that I originally developed for the purposes of running statistical tests over large GWAS cohorts in R: RegParallel.
Another program that can run a test directly on your VCF data is SnpSift CaseControl.
For a Manhattan plot, you will require:
You then just need the qqman package:
library(qqman)
manhattan(
subset(temp, select=c(SNP, CHR, BP, P)),
main="Allelic model",
suggestiveline=-log10(0.0001),
genomewideline=-log10(5e-08),
col=c("springgreen4", "firebrick"),
chrlabs=c(1:22, "X", "Y", "MT"),
ylim=c(0,15))
legend("topright", cex=0.8, title="Significances", c("P<0.0001", "FDR (P<5.2E-08)"), fill=c("blue", "red"))
Hello!
Thank you for your response.
I have 53 VCF files of patients (53 files for 53 patients)
I checked the A: SNP dataset and Z Score, I don't know how to find out allele frequency.
My dataset consists of Chromosome position, RSID, ref genome, mutated genome, gene name , amino acid change, 1000 genome exome aggregation database for all,ExAC ALL.
If you just have 53 patients, which statistical test do you want to perform? ...or is it that you just want to calculate the allele frequency?
Your 'dataset' is all of the data from the patients combined together? Can you paste a few rows from the data?
I have been given 53 patients now to practice and learn. I want to know about PLINK or R whichever is easier to work with.
No,, it is not just the allele frequency that I want to calculate, ultimately I have to plot the Manhattan plot.
I have separate vcf files for each of 53 patients. I also have a combined excel file for all patients together.
ID SID CH POS RSid Ref Alt GeneName AAChange 1000g2015aug_all ExAC_ALL
1 IS1 1 13656 . CAG C DDX11L1 . . .
2 IS1 1 761957 rs59038458 A AT LINC00115 . 0.265176 .
3 IS1 1 874950 rs6143081 T TCCCTGGAGGACC SAMD11 . 0.729832 .
this is of 1 patient only which is identified as IS1
PS: I cannot share actual data.
Thanks.
Regarding the Manhattan plot, this plot is used to show the distribution of p-values, genome-wide. If your dataset is just 53 patients with no controls, then you cannot derive any p-values. So... who told you to generate a Manhattan plot? Do they understand that your dataset is 100% patients? Perhaps they meant something else.
The table that you have looks like output from the common program called 'ANNOVAR'. How did you obtain this table? The person who gave it to you obviously has done some pre-processing of the data.
My first step, if I were you, would be to merge all individual VCFs into a single VCF file. You can do this with bcftools merge --merge none [VCF] [VCF] ... [VCF]
"If your dataset is just 53 patients with no controls, then you cannot derive any p-values. So... who told you to generate a Manhattan plot? Do they understand that your dataset is 100% patients? Perhaps they meant something else" That's why I mentioned that I need to compare my data with the healthy controls and how do I do that!
"*The table that you have looks like output from the common program called 'ANNOVAR'. How did you obtain this table? *" Yes, it is a preprocessed data (an Excel file), However I also have vcf file. I do not how which tool is appropriate to open it with so that I can read the data.
"merge all individual VCFs into a single VCF file. You can do this with bcftools merge --merge none [VCF] [VCF] ... [VCF]" Thanks! Also,Do I need to install bcftools only, to run this from cmd ? can I work it on windows 10 or Ubuntu should be preferably better?
After I merge all the Vcf files into 1 VCF file, what is the next step?
Thank you for patiently answering :)
That's why I mentioned that I need to compare my data with the healthy controls and how do I do that!
Cool / vale, but, which controls? You have none. Do you mean, for example, 1000 Genomes Controls?
Yes, it is a preprocessed data (an Excel file), However I also have vcf file. I do not how which tool is appropriate to open it with so that I can read the data.
So, you already have a single VCF file that contains all samples?
Thanks! Also,Do I need to install bcftools only, to run this from cmd ? can I work it on windows 10 or Ubuntu should be preferably better?
Ubuntu will be better.
Is your supervisor not a bioinformatician?
If you want to obtain all 1000 Genomes data in VCF format, then you can follow the first few steps, here: Produce PCA bi-plot for 1000 Genomes Phase III - Version 2
Once you have both datasets in separate VCFs, that would be a great first step.
Note that the 1000 Genomes data is ~15 gigabytes
Sadly No. My entire lab is a Bionano research area where my colleagues and my supervisors specialize in proteomics and immunology and bionanoscience. My postgraduate thesis is this and I can use any help I can get. Thanks!
I will get back to the thread when I am done with the merging up vcf files. Thank you again!
Hi
I merged the 53 vcf files into one using the command: bcftools merge --merge none [VCF] [VCF] ... [VCF]
however my concerns are: a) is the command same for even 200 or 1000 samples? Or, is there a way to merge all the 1000 vcf files into one without typing all the 1000 file names into the command?
b) This command merges all the data from all the vcf files mentioned in the command. However, is there a command which calls in specific columns and then merges those particular columns only for all the n number of vcf files?
Thank you!
Ye, you can merge any number of VCFs using BCFtools merge. What you can do is do it programmatically, though.
If you have your VCF files in a single column in VCF.list:
command="bcftools merge -Ob -m none " ;
while read VCF
do
command="${command}"" ""${VCF}" ;
done < VCF.list
command="${command}"" -o ProjectMerge.bcf" ;
echo `$command` ;
bcftools index ProjectMerge.bcf ;
This will merge any number of VCFs, i.e., all of the VCFs listed in VCF.list. Note that, in the penultimate echo command, there are back-ticks used, not apostrophes.
b) This command merges all the data from all the vcf files mentioned in the command. However, is there a command which calls in specific columns and then merges those particular columns only for all the n number of vcf files?
I am not sure what you mean?
b) This command merges all the data from all the vcf files mentioned in the command. However, is there a command which calls in specific columns and then merges those particular columns only for all the n number of vcf files?
the vcf files contain many columns, from which i only want to filter out these columns: ID | SID |CH | POS | RSid | Ref | Alt | GeneName | AAChange | 1000g2015aug_all | ExAC_ALL
command="bcftools merge -Ob -m none " ;
while read VCF
do
command="${command}"" ""${VCF}" ;
done < VCF.list
command="${command}"" -o ProjectMerge.bcf" ;
echo `$command` ;
bcftools index ProjectMerge.bcf ;
Do you mind explaining this code to me? I would like to understand and not blindly follow. Also, Do I write this code on terminal? Or is it a python code?
If your VCF files has columns called ID, SID, CH, POS, et cetera, then it is not a VCF file. Please review the VCF specification: http://www.internationalgenome.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41/
Maybe you just mean that you want to extract these from the INFO column in the VCF?
Take a look at your VCF header:
bcftools view -h Greek.bcf
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##pipeline_version=DVDS_v1.4.7.5
##resource_version=v1.0.13
##FILTER=<ID=INDEL_SPECIFIC_FILTERS,Description="QD < 2.0 || ReadPosRankSum < -20.0 || InbreedingCoeff < -0.8 || FS > 200.0">
##FILTER=<ID=LowQual,Description="Low quality">
##FILTER=<ID=VQSRTrancheSNP99.00to99.90,Description="Truth sensitivity tranche level for SNP model at VQS Lod: -4.0314 <= x < -0.6424">
##FILTER=<ID=VQSRTrancheSNP99.90to100.00+,Description="Truth sensitivity tranche level for SNP model at VQS Lod < -2467.098">
##FILTER=<ID=VQSRTrancheSNP99.90to100.00,Description="Truth sensitivity tranche level for SNP model at VQS Lod: -2467.098 <= x < -4.0314">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP Membership">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
##INFO=<ID=Dels,Number=1,Type=Float,Description="Fraction of Reads Containing Spanning Deletions">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
##INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with at most two segregating haplotypes">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
##INFO=<ID=MQ0,Number=1,Type=Integer,Description="Total Mapping Quality Zero Reads">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=NEGATIVE_TRAIN_SITE,Number=0,Type=Flag,Description="This variant was used to build the negative training set of bad variants">
##INFO=<ID=POSITIVE_TRAIN_SITE,Number=0,Type=Flag,Description="This variant was used to build the positive training set of good variants">
All of those key-value pairs can be extracted with a single command:
bcftools query -f '%FS\t%AC\t%DP\t%HaplotypeScore\n' Greek.bcf | head
0 2 5012 0.1734
0 0 3234 0.4606
0 0 3271 0.549
0 7 26770 0
7.384 0 8629 0.7159
1.182 0 7623 1.9923
0 0 1767 0.4625
0 32 703 0
0 14 890 0
3.854 3 870 0.035
Regarding the code that I shared: It is essentially constructing the bcftools merge command as a string variable. In each loop, a new file is appended to the string. By the final loop, the string is just:
"bcftools merge -Ob -m none VCF VCF VCF VCF VCF ... VCF -o ProjectMerge.bcf"
This is then executed by the echo command using back-ticks `
Hi Kevin,
Thank you so much for the help. I understood now how to read the vcf files and how to merge multiple files with a single command into one.
So now, I download the GrCH38 fastq files from internationalgenome of the controls. Do i convert them to bcf files like you mentioned in your command here:
for chr in {1..22}; do
bcftools norm -m-any --check-ref w -f human_g1k_v37.fasta \
ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | \
bcftools annotate -x ID -I +'%CHROM:%POS:%REF:%ALT' |
bcftools norm -Ob --rm-dup both \
> ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.bcf ;
bcftools index ALL.chr"${chr}".phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.bcf ;
done
if yes, why are we converting them to bcf files, and not vcf or other format? I know bcf is binary format, but i wanna know the significance.
Hey, for these 'smaller' questions, like differences between VCF and BCF, you may want to search for the information in a search engine.
Relating to the code (above), that code normalises the VCF files and then converts them to BCF. BCF is not BAM format. By 'normalise', we mean that it left-aligns indels according to the reference genome (human_g1k_v37.fasta) and also splits multi-allelic variants to multiple rows / lines. If you have data in GRCh38 co-ordinates, then you cannot use human_g1k_v37.fasta because human_g1k_v37.fasta is GRCh37.
If you want to use the code above, you can delete the --check-ref w -f human_g1k_v37.fasta part.
So now, I download the GrCH38 fastq files from internationalgenome of the controls. Do i convert them to bam files like you mentioned in your command here:
From where have you downloaded the files?
It was an error in writing bam, I meant bcf.
My question was simply that if I have the fasta file of the controls, why am I not converting them to vcf file format? I know bcf is the binary format of vcf and is efficient to process. So, then do I convert the vcf file for cases to bcf as well?
PS: I am a Student learning this from scratch. I don't think I should shy away from asking any question smaller or bigger.
Ahora será mejor abrir una nueva pregunta. En la nueva pregunta, relata sobre lo que ya has hecho y lo que quieres hacer con los archivos // dados que tiene.
By Not answering my question you just pushed me to do better and take this as a challenge. The world doesn't stop when one person who knows the subject doesn't wanna help others learn because they ask smaller questions which can be easily googled. Sometimes, an explanation from a teacher (in this case, you) is far better understandable and helpful rather than reading those from some papers or books online. That's why they have teachers and professors to explain those writings from books and search engines to students who want to learn. I hope you too learn how to be patient with students who have absolutely no idea about a subject but want to learn.
At this point, I don't want your help anymore.
Thank you for the arrogance!
Our work on this website is 100% voluntary. We receive no payment for coming here and we give time from our days to help users, most of whom are new to the field of bioinformatics. Virtually all of us have full-time jobs and families.
Each day, there are multiple new posts and many new questions asked. I have already answered many of your questions for 10 days in this thread. It is your responsibility as a student to take on board my help and then learn from it.
Mañana será otro día.
Te deseo lo mejor / Saludos.
It was your own incorrect interpretation that I was thinking along those lines.
By Not answering my question you just pushed me to do better and take this as a challenge.
Then, your experience here has been positive and that is what I was hoping that you would do.
Everyone who participates on this website by answering questions are doing this voluntarily. We do not get paid or benefit from this. We all have other responsibilities in our professional and personal lives. We are not here to act as your professor or full-time teacher.
I understand that learning bioinformatics can be challenging and frustrating at times. But please do not direct your frustration towards people who are just trying to balance their own responsibilities with contributing to the bioinformatics community.
Please show the preceeding commands that you have used.
Thank you!
고맙습니다
¡Gracias!
Log in to answer this question.
"If your dataset is just 53 patients with no controls, then you cannot derive any p-values. So... who told you to generate a Manhattan plot? Do they understand that your dataset is 100% patients? Perhaps they meant something else" That's why I mentioned that I need to compare my data with the healthy controls and how do I do that!
"*The table that you have looks like output from the common program called 'ANNOVAR'. How did you obtain this table? *" Yes, it is a preprocessed data (an Excel file), However I also have vcf file. I do not how which tool is appropriate to open it with so that I can read the data.
"merge all individual VCFs into a single VCF file. You can do this with bcftools merge --merge none [VCF] [VCF] ... [VCF]" Thanks! Also, how? I mean to ask, is this cmd? Do i need special tool to run this command?
Thank you for patiently answering :)
Please use
ADD COMMENT/ADD REPLYwhen responding to existing posts to keep threads logically organized.Re-post this comment at its logical location in @Kevin's answer below and then come back to delete this. You can use " icon in the edit tool bar to quote someone else' text.
Hello jaybee!
We believe that this post does not fit the main topic of this site.
If you still require help with this, then please open a new question that is very specific
For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.
If you disagree please tell us why in a reply below, we'll be happy to talk about it.
Cheers!