Question: (Closed) New to Whole Exome Sequencing and Manhattan Plot and I need to understand how to compare patients' data with healthy control and find p value and then plotting Manhattan plot.
0
gravatar for jaybee
15 days ago by
jaybee10
South Korea
jaybee10 wrote:

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.

ADD COMMENTlink modified 14 days ago • written 15 days ago by jaybee10

"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 :)

ADD REPLYlink written 14 days ago by jaybee10

Please use ADD COMMENT/ADD REPLY when 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.

ADD REPLYlink written 14 days ago by genomax62k

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!

ADD REPLYlink written 1 day ago by Kevin Blighe37k
0
gravatar for Kevin Blighe
14 days ago by
Kevin Blighe37k
Republic of Ireland
Kevin Blighe37k wrote:

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:

  • ID (e.g. rs ID)
  • CHR (chromosome)
  • BP (base position)
  • P (p-value)

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"))

g

ADD COMMENTlink written 14 days ago by Kevin Blighe37k

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.

ADD REPLYlink written 14 days ago by jaybee10

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?

ADD REPLYlink written 14 days ago by Kevin Blighe37k

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.

ADD REPLYlink modified 2 days ago by RamRS20k • written 14 days ago by jaybee10

https://photos.app.goo.gl/Vt4Lyxj7jQzES1hw6

ADD REPLYlink written 14 days ago by jaybee10
1

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]

ADD REPLYlink written 14 days ago by Kevin Blighe37k

"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 :)

ADD REPLYlink modified 14 days ago • written 14 days ago by jaybee10
1

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.

ADD REPLYlink written 14 days ago by Kevin Blighe37k

Yes, healthy controls from 1000 genome.org

after I merge all the vcf files into 1 vcf file, what is my next step?

Thanks!

ADD REPLYlink written 13 days ago by jaybee10

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

ADD REPLYlink written 13 days ago by Kevin Blighe37k
1

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!

ADD REPLYlink written 13 days ago by jaybee10

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!

ADD REPLYlink written 8 days ago by jaybee10

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?

ADD REPLYlink modified 8 days ago • written 8 days ago by Kevin Blighe37k

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?

ADD REPLYlink modified 7 days ago • written 7 days ago by jaybee10
1

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 `

ADD REPLYlink modified 7 days ago • written 7 days ago by Kevin Blighe37k

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.

ADD REPLYlink modified 6 days ago • written 6 days ago by jaybee10

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.

ADD REPLYlink written 6 days ago by Kevin Blighe37k

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?

ADD REPLYlink written 6 days ago by Kevin Blighe37k

http://www.internationalgenome.org/

ADD REPLYlink written 6 days ago by jaybee10

Also I know the difference between VCF and BCF and I know BCF is not BAM. Also, I know GRCh 38 is human_g1k_v38.fasta and the reference is hg38

ADD REPLYlink written 6 days ago by jaybee10

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.

ADD REPLYlink modified 6 days ago • written 6 days ago by jaybee10
1

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.

ADD REPLYlink written 3 days ago by Kevin Blighe37k

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!

ADD REPLYlink modified 3 days ago • written 3 days ago by jaybee10
2

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.

ADD REPLYlink written 3 days ago by Kevin Blighe37k
2

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.

ADD REPLYlink written 3 days ago by Damian Kao15k

If you think we do not benefit from this, then I really don't have anything to say. We all benefit from learning and/or teaching.

ADD REPLYlink written 1 day ago by jaybee10

Also, I do not find learning bioinformatics frustrating. That is a very poor usage of the word "frustrating". Please do not use such words and try and change what I wanted to say into your own perspective.

Thank you!

ADD REPLYlink written 1 day ago by jaybee10

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.

ADD REPLYlink modified 3 days ago • written 3 days ago by Kevin Blighe37k

And, I did. But nothing gives you any right to tell anyone if their question is small or big! :)

ADD REPLYlink written 1 day ago by jaybee10

It was your own incorrect interpretation that I was thinking along those lines.

ADD REPLYlink modified 1 day ago • written 1 day ago by Kevin Blighe37k

after I run

echo `$command` ;

It gives me error message:

Failed to open vcf.list: unknown file type

ADD REPLYlink written 1 day ago by jaybee10

Please show the preceeding commands that you have used.

Thank you!

고맙습니다

¡Gracias!

ADD REPLYlink modified 1 day ago • written 1 day ago by Kevin Blighe37k
Please log in to add an answer.
The thread is closed. No new answers may be added.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2010 users visited in the last hour