AlanFung:LabNotes/ASE/2010/2012-7-18
Jump to navigation
Jump to search
Capturing Protocol[edit]
[[1]]
Analysis of Gimelbrant Lab Samples[edit]
- Data Stored in seqStore/120712_HL125, transfer my data to Gimelbrant/ase_06262012/tmp directory
- Transfer to my own directory
- Create folder gimelbrant/ase_06262012
- Make a directory for each samples (bash mkdir.txt)
Combine files and transfer to individual folders[edit]
- Run combine.txt to combine data and separate them into individual folders
- Run move.txt to transfer files into individual folders
AseBowtieSam[edit]
- Run aseBowtieSam.pl
- chmod +x../aseBowtie.pl to make is executable
- make sure the directories are up to date by running each line manually if not up to date edit the script
my $mapper = "/home/kunzhang/softwares/bowtie-0.11.3/bowtie"; * change to my $mapper = "/home/kunzhang/softwares/bowtie-0.12.8/bowtie";
my $samtool = "/home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools";
my $readTrimmer="/home/kunzhang/Exome/Solexa/scripts/trimTerminalMismatches.pl"; *Cannot find this one *change to my $readTrimmer="/home/kunzhang/Exome/Data/scripts/trimTerminalMismatches.pl";
my $bowtie2sam = "/home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl";
my $snpFilter = "/home/kunzhang/MIP/pileup2basecalls.pl";
my $in_file = $ARGV[0];
my $template = "/home/kunzhang/HsGenome/hg18.masked/hg18.masked"; my $fai_file = "/home/kunzhang/HsGenome/hg18.masked/HsGenome18.masked.fai";
- Run the script with the .txt file
- after editing the script, copy script to cDNA-DF2
- run script
- Return with error
/usr/local/bin/maq sol2sanger cDNA-DF2.txt cDNA-DF2.fastq Can't exec "/usr/local/bin/maq": No such file or directory at aseBowtieSam.pl line 28. Warning: Could not open read file "cDNA-DF2.fastq" for reading; skipping... sh: cannot open cDNA-DF2.bowtie.out: No such file Can't open cDNA-DF2.bowtie.out.trimmed: No such file or directory at /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl line 19. Can't use an undefined value as an ARRAY reference at /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl line 48.
- It appears that maq is uninstalled, Athurva will help me install the software
Shell to run script on all samples[edit]
*for f in gDNA-H14 ChIPK36me3-P10 ChIPK36me3-DF2 ChIPK27me3-P10 ChIPK27me3-DF2 cDNA-P10E3 cDNA-DF2 gDNA-GM gDNA-P10E3 gDNA-DF2 gDNA-H0 gDNA-GmgDNA:H0gDNA::1:1 do cd $f ../aseBowtieSam.pl $f.txt >$f.log cd .. done
- Run is complete
- What this script does? Looking at the output it generates
.bowtie.basecalls.txt .pileup .sorted.bam
- Going back to Athurava's lecture notes
- Pileup is the consensus sequence, it describes the base-pair information at each chromosomal position
Pileup File Columns 1. Chromosome 2. Position 3. reference base 4. Call 5. Consensus quality
P(call column is incorrect) Coverage sum of quality score at that position
Analysis for each sample[edit]
Total reads[edit]
- wc -l .txt >> .txt
- go into the fastq .txt files and divide the total lines by 4
- run it using script
- copy to excel data to text column /4
Mapping Rate[edit]
- /../home/kunzhang/softwares/samtools-latest/samtools flagstat cDNA-DF2.sorted.bam
- line 4 is the mapped reads
- sed -n 4p $f.log2 >> ../mapreadcount.txt
- to print line 4 into text file
#SNPs called[edit]
- wc -l $f.bowtie.basecalls.txt >> ../snpcount.txt
- Basically count the lines in the basecalls.txt file
#Het called[edit]
- Kun told me i can use the IF function in excel and count the heterozygous SNP
- Used IF function for column C , =IF(LEFT(C1)=RIGHT(C1),0,1)
- If the two base are the same returns 0
- If the two bases are different returns 1
- Get the total of the column and you get the total number of heterozygous SNP
Data Summary[edit]
Sample | Total Number of Reads | Mapping Reads | Mapping Rate | #SNPs called | #Het | % Het |
gDNA-H7 | 7,214,167 | 4,946,303 | 68.56 | 28,280 | 2954 | 10.45 |
gDNA-H14 | 9,686,002 | 8,221,702 | 84.88 | 20,536 | 2502 | 12.18 |
ChIPK36me3-P10 | 8,811,298 | 7,223,409 | 81.98 | 25,452 | 2451 | 9.63 |
ChIPK36me3-DF2 | 9,243,720 | 7,649,171 | 82.75 | 26,507 | 2660 | 10.04 |
ChIPK27me3-P10 | 9,149,672 | 7,490,869 | 81.87 | 25,320 | 2475 | 9.77 |
ChIPK27me3-DF2 | 8,125,103 | 6,607,398 | 81.32 | 26,137 | 2687 | 10.28 |
cDNA-P10E3 | 6,012,037 | 5,207,362 | 86.62 | 19,096 | 2451 | 12.84 |
cDNA-DF2 | 4,306,791 | 3,756,778 | 87.23 | 18,854 | 2405 | 12.76 |
gDNA-GM | 6,920,955 | 5,504,556 | 79.53 | 30,270 | 4188 | 13.84 |
gDNA-P10E3 | 7,540,086 | 6,126,998 | 81.26 | 30,887 | 4223 | 13.67 |
gDNA-DF2 | 6,836,212 | 5,328,454 | 77.94 | 30,249 | 4148 | 13.71 |
gDNA-H0 | 8,268,386 | 6,506,594 | 78.69 | 30,809 | 3590 | 11.65 |
gDNA-GmgDNA:H0gDNA::1:1 | 9,113,177 | 7,336,107 | 80.50 | 29,753 | 4065 | 13.66 |