1.PE151测序需要多少M PE Read-peir 覆盖？一个10Gbp大小的genome 99%以上的区域达到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基因的变异，说明变异位点是的数目（纯和和杂合）。
#提取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 软件 ，设置 -K参数为35对该数据 进行de novo组装 ，并 画出组装结果序列从长到短的度累积曲线图。
计算组装结果的 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是蛋白质TP53在19种动物中的同源序列集：文件program/mafft是一款易用的多序列比对软件；文件program/FastTreeMP是一款快速构建系统发育数的软件;
构建系统发育树。
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()
image.png
7.文件program/trf是 Tandem Repeat Finder软件,利用Tandem Repeat Finder软件分析data/samll_genome.fasta.gz中可能的端粒序列，回答该基因组的端粒序列重复单元为？
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
image.png
端粒是短的多重复的非转录序列（TTAGGG）及一些结合蛋白组成特殊结构
端粒序列:CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAA 重复单元：CCTAAA对应互补链TTTAGGG 人的端粒是短的多重复的非转录序列（TTAGGG）及一些结合蛋白组成特殊结构。端粒是存在于真核细胞线状染色体末端的一小段DNA-蛋白质复合体,端粒通常由富含鸟嘌呤核苷酸（G）的短的串联重复序列组成，端粒DNA是由简单的DNA高度重复序列组成的，染色体末端沿着5'到3' 方向的链富含GT。
8.data/small_gene_set.gff.gz是某植物基因信息文件，data/genome.fasta.gz是该植物的基因组序列文件，而program/gffread是gff注释信息转换软件，-E校验gff文件并转换为标准格式；-x参数参数提取基因的编码序列；Program/cds2aa.pl文件是一个Perl小程序；
提取CG07320.UN 的注释信息保存在666.gff中，利用gffread软件提取small_genoeme.fasta CDS序列翻译出对应氨基酸序列，基因最长转录物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
image.png
[kmer_Freq]$ cat read_list
newBGIseq500_1.fq.gz
newBGIseq500_2.fq.gz
1、accession号：XP_008224950.1。物种： Prunus mume
2、该相近蛋白具有一个结构域，名称为WD40 super family，基因名为Wrap53，主要的molecular funtion为RNA binding（GO:0003723）和protein binding（GO:0005515）
3、搜索telomerase Cajal body protein 1 isoform X2发现0篇相关文献
搜索 telomerase Cajal body protein 1 发现21篇相关文献
1、最可能是Oryza punctata 斑点野生稻;
2、相差9个位点。只有695位点的G->T为最可能的真实差异位点。464，514，516位点为杂合峰，695为单峰，700位以后的QV值都非常低，说明测序质量低，不是真实的突变位点。
操作过程进入NCBI网站，选择tools，选择Basic local alignment search tools（BLAST）,然后选择nucleotide blast。然后点击浏览上传rbcL.seq。输入job title。点击BLAST提交任务。Sequences producing significant alignments:的第一条为比分得分最高的，点击Aaccession处对于的id。查看详细比对情况和sanger测序峰图。
网友评论