1.PE151測序需要多少M PE Read-peir 覆蓋？一個10Gbp大小的genome 99%以上的區(qū)域達到10x也糊。
# 10*10^9 *0.99*10 = 150 *2*x
需要330M
R代碼如下
l<-150*2
total<-10*1000000000*0.99*10
reads<-total/l
final<-reads/1000000
2.人類基因組計劃和后基因組計劃
人類基因組計劃圖譜:
遺傳圖譜（genetic map）,物理圖譜（physical map）,序列圖譜,基因圖譜
后基因組計劃:
HapMap:國際人類基因組單倍型圖計劃
旨在識別人類基因組中所有功能因素的DNA百科全書（簡稱ENCODE）
1000 genome project:千人基因組計劃
ICGC:癌癥基因組計劃
HEP:Human Epigenome Project,人類表觀基因組計劃
3.計算vcf文件中Qual值分別情況，畫圖技矮；
提取TP53基因的變異吟策，說明變異位點是的數(shù)目（純和和雜合）。
#提取Qual值
1.vcftools --gzvcf chr17.vcf.gz --site-quality --out Q
or
2. less -S chr17.vcf.gz | grep '^chr17' | awk -F "\t" '{print $2"\t"$6}' > chr17.qual
#畫圖
1.vcftools --gzvcf chr17.vcf.gz --chr chr17 --from-bp 7668402 --to-bp 7687550 --het -c | head
or
2. zcat chr17.vcf.gz | grep '^chr17' | awk -F "\t" '{if ($2>=7668402 && $2<=7687550){print $0}}' | wc -l
INDV O(HOM) E(HOM) N_SITES F
27DMBDM4YT 24 16.0 24 1.00000
7XKZJA3JWX 0 16.0 24 -2.00000
APRDKV0LDS 24 16.0 24 1.00000
awk -F "\t" '{print $11}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if ($1==$2){print $0}}'
完整代碼
$$
## qual.R
png("QualStatistics.png")
a <- read.table("Q.lqual",header = TRUE)
q <- a[,3]
hist(q,main = "Qual Hist")
dev.off()
vcftools --gzvcf chr17.vcf.gz --site-quality --out Q
Rscript qual.R
#vcftools --gzvcf chr17.vcf.gz --chr chr17 --from-bp 7668402 --to-bp 7687550 --het -c | head
#vcftools --gzvcf chr17.vcf.gz --bed tp54.bed --recode --out tp54
vcftools --gzvcf chr17.vcf.gz --chr chr17 --from-bp 7668402 --to-bp 7687550 --recode --out tp54
less -S tp54.recode.vcf | grep '^chr17' | wc -l
for i in {10..12};do
cat tp54.recode.vcf | grep '#CHROM' |awk -F "\t" '{print $'$i'}'
echo -e "hom\t `cat tp54.recode.vcf | grep '^chr17' |awk -F "\t" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1==$2){print $0}}' | wc -l`"
echo -e "het\t `cat tp54.recode.vcf | grep '^chr17' |awk -F "\t" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1!=$2){print $0}}' | wc -l`"
done
4.下載并安裝 SOAPdenovo 軟件 件相，設(shè)置 -K參數(shù)為35對該數(shù)據(jù) 進行de novo組裝 ，并 畫出組裝結(jié)果序列從長到短的度累積曲線圖氧苍。
計算組裝結(jié)果的 N50
下載地址：https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/bin/r240/SOAPdenovo2-bin-LINUX-generic-r240.tgz/download
tar -zxf SOAPdenovo2-bin-LINUX-generic-r240.tgz
or
ftp://public.genomics.org.cn/BGI/SOAPdenovo2/Assemblathon1_pipeline.tgz
#從大到小累計圖
png("length.png")
data<-read.table("zzz",head=F)
colnames(data)<-c("length")
pdata<-table(data$length)
pframe<-data.frame(as.numeric(names(pdata)),as.numeric(pdata))
colnames(pframe)<-c("length","num")
rankData<-pframe[order(pframe[,1],decreasing=T),]
mulData<-data.frame()
sum=0
for(i in 1:nrow(rankData)){
sum=sum+rankData[i,2]
row=cbind(rankData$length[i],sum)
mulData=rbind(mulData,row)
}
colnames(mulData)<-c("length","num")
plot(mulData$length,mulData$num/sum,main="序列長度累積曲線圖",xlab="長度",ylab="累計比例",type="l")
axis(side=2,seq(from=0,to=1,by=0.1),)
#axis(side=1,seq(from=36000,to=300000,by=10000),)
dev.off()
# cat config
max_rd_len=100
[LIB]
avg_ins=450
asm_flag=3
map_len=32
pair_num_cufoff=3
reverse_seq=0
rank=1
q1=newBGIseq500_1.fq.gz
q2=newBGIseq500_2.fq.gz
#cat work.sh
SOAPdenovo-63mer_v2.0 all -s config -K 35 -D 1 -o species > log 2> err
perl 05_result.pl species.scafSeq
Rscript plot.R
5.TP53.pep.fa是蛋白質(zhì)TP53在19種動物中的同源序列集：文件program/mafft是一款易用的多序列比對軟件夜矗；文件program/FastTreeMP是一款快速構(gòu)建系統(tǒng)發(fā)育數(shù)的軟件;
構(gòu)建系統(tǒng)發(fā)育樹。
mafft TP53.pep.fa > TP53.fa
FastTreeMP TP53.fa >TP53.newick
Rscript TP53_phylogeny.R
cat TP53_phylogeny.R
png("TP53_phylogeny.png")
#library("ggplot2")
library("ape")
tree <-read.tree(file="TP53.newick")
plot(tree)
dev.off()
#可視化：
1,生成的Newick tree格式文件可視化工具：http://www.trex.uqam.ca/index.php?action=newick
2,可視化工具可以用PHYLIP programs（下載地址
http://evolution.genetics.washington.edu/phylip/getme.html)
3,還有R中的library(ape)
6.畫kmer分布圖
png("kmer_distibution.png")
a <- read.delim("kmer_freq_distribution",head=F)
l<-seq(1,nrow(a),by=1 )
plot(x=l,y=a[,2],type ="l",col ="red",lwd=2,xlim=c(0,200),ylim=c(0,1.5),main="kmer頻率分布圖",xlab="depth",ylab="Frequency(%)")
#text(40,0.05,"38")
dev.off()
7.文件program/trf是 Tandem Repeat Finder軟件,利用Tandem Repeat Finder軟件分析data/samll_genome.fasta.gz中可能的端粒序列候引，回答該基因組的端粒序列重復(fù)單元為侯养？
gunzip small_genome.fasta.gz
example/program/trf small_genome.fasta 2 7 7 80 10 50 2000 -d -h
less -S small_genome.fasta.2.7.7.80.10.50.2000.dat
端粒是短的多重復(fù)的非轉(zhuǎn)錄序列（TTAGGG）及一些結(jié)合蛋白組成特殊結(jié)構(gòu)
端粒序列:CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAA 重復(fù)單元：CCTAAA對應(yīng)互補鏈TTTAGGG 人的端粒是短的多重復(fù)的非轉(zhuǎn)錄序列（TTAGGG）及一些結(jié)合蛋白組成特殊結(jié)構(gòu)。端粒是存在于真核細胞線狀染色體末端的一小段DNA-蛋白質(zhì)復(fù)合體,端粒通常由富含鳥嘌呤核苷酸（G）的短的串聯(lián)重復(fù)序列組成澄干，端粒DNA是由簡單的DNA高度重復(fù)序列組成的逛揩，染色體末端沿著5'到3' 方向的鏈富含GT。
8.data/small_gene_set.gff.gz是某植物基因信息文件麸俘，data/genome.fasta.gz是該植物的基因組序列文件辩稽，而program/gffread是gff注釋信息轉(zhuǎn)換軟件，-E校驗gff文件并轉(zhuǎn)換為標(biāo)準(zhǔn)格式从媚；-x參數(shù)參數(shù)提取基因的編碼序列逞泄；Program/cds2aa.pl文件是一個Perl小程序；
提取CG07320.UN 的注釋信息保存在666.gff中拜效，利用gffread軟件提取small_genoeme.fasta CDS序列翻譯出對應(yīng)氨基酸序列喷众，基因最長轉(zhuǎn)錄物mRNA等3個密碼子及氨基酸名稱。
zcat small_gene_set.gff.gz | grep 'CG07320' > 666.gff
gunzip small_genome.fasta.gz
program/gffread 666.gff -E -o- >666.e.gff
program/gffread 666.e.gff -g small_genome.fasta -x 666.cds.fa -o out.gff
perl ../../program/cds2aa.pl 666.cds.fa >666.pep.fa
less -S 666.pep.fa
[kmer_Freq]$ cat read_list
newBGIseq500_1.fq.gz
newBGIseq500_2.fq.gz
1紧憾、accession號：XP_008224950.1到千。物種： Prunus mume
2、該相近蛋白具有一個結(jié)構(gòu)域赴穗，名稱為WD40 super family憔四，基因名為Wrap53膀息，主要的molecular funtion為RNA binding（GO:0003723）和protein binding（GO:0005515）
3、搜索telomerase Cajal body protein 1 isoform X2發(fā)現(xiàn)0篇相關(guān)文獻
搜索 telomerase Cajal body protein 1 發(fā)現(xiàn)21篇相關(guān)文獻
1了赵、最可能是Oryza punctata 斑點野生稻;
2潜支、相差9個位點。只有695位點的G->T為最可能的真實差異位點柿汛。464冗酿，514，516位點為雜合峰苛茂，695為單峰已烤，700位以后的QV值都非常低，說明測序質(zhì)量低妓羊，不是真實的突變位點。
操作過程進入NCBI網(wǎng)站稍计，選擇tools躁绸，選擇Basic local alignment search tools（BLAST）,然后選擇nucleotide blast。然后點擊瀏覽上傳rbcL.seq臣嚣。輸入job title净刮。點擊BLAST提交任務(wù)。Sequences producing significant alignments:的第一條為比分得分最高的硅则，點擊Aaccession處對于的id淹父。查看詳細比對情況和sanger測序峰圖。
