- 使用misa.pl 查找ssr位点
生成的misa文件有
wc -l tsh_wgs.fna.misa
127131 tsh_wgs.fna.misa
#生成的misa文件行数,但statistics文件里显示Total number of identified SSRs: 146445
Number of SSRs present in compound formation: 19315
即总共的SSRs 146445 = 127131(行数)- 1(表头)+19315(复合卫星数量)=146445 也即复合卫星计算为2个SSR
3. p3_in.pl 将misa.pl生成的misa文件处理成primer3 可以识别的内容
网上资料说需要改p3_in.pl 里的内容,将print OUT 后内容更改成primer可以识别的内容,之后运行p3_in.pl ,生成的p3in文件内容为空。之后发现可能是因为基因组序列文件> 后的序列名称含有空格,而misa.pl程序会用下划线补全空格。可能是因为无法对应misa文件与基因组的序列名称,导致生成空文件。
4. 将基因组序列名称只保留是scaffold内容,序列名称无空格。
sed -i 's/>.Taishanhong />/;s/, .//' fasta文件名称,将Taishanhong以前的内容替换为空,将 ","以后的内容替换为空
5. 可以运行 p3_in.pl 不过生成的 p3in文件有204GB,好奇怪。也没搞懂 p3in文件里的SEQUENCE_TEMPLATE=是如何得到的。
6. 查看生成的p3in文件,p3_in.pl 的print OUT 内容只有这四行,
其中
6.1 SEQUENCE_ID=scaffold1_1 表示scaffold1 上查找到的第一个位点,
想看一下scaffold1总共查到到多少个位点,使用
grep "scaffold1_" tsh_wgs.fna.p3in |wc -l
grep -c "scaffold1_" tsh_wgs.fna.p3in
感觉很慢,没有出现结果,直接 grep "scaffold1_" tsh_wgs.fna.p3in 查看scaffold上有3793个位点。
6.2 SEQUENCE_TEMPLATE= 猜测模板序列可能是scaffold1 整条序列的碱基。
sed -n 2p tsh_wgs.fnarow1-4.p3in |awk -F "=" '{print $2}' | wc -m
6916724 #p3inp3in文件里 SEQUENCE_ID=scaffold1_1 的模板序列碱基数量
需要统计一下基因组数据里scaffold1的碱基数量,fna文件里scaffold的碱基是每行的数量固定,有很多行,就需要提取scaffold1 与 下一个scaffold之间的内容并计算字符数。
grep -n "scaffold" tsh_wgs.fna |less -SN
1 1:>scaffold1
2 86462:>scaffold10 #即scaffold1 的内容在第1行于第86462之间,也即从第2行到 86461行。tsh_wgs.fna 文件中scaffold1共有86461-2+1= 8640 行。
3 134759:>scaffold100
sed -n '2,86461p' tsh_wgs.fna | wc -l
86460 #提取了86460行碱基
sed -n '2,86461p' tsh_wgs.fna | wc -m
7003183 #与p3in 的碱基数量不一致,可能wc -m计算的不仅仅是碱基字母数量
sed -n '2p' tsh_wgs.fna | wc -m
81 #即第二行有81个字符,但是tsh_wgs.fna文件里第二行只有80个碱基,
sed -n '2p' tsh_wgs.fna | awk -F "" '{print $80}'
G #为第二行的最后一个碱基
了解到第81个字符为行尾符。
说明上边的计算没错。tsh_wgs.fna scaffold1里的
字符数 7003183 - 总共的行尾符数量 86460(即总行数86460) = 总碱基数
=p3in文件里的 SEQUENCE_ID=scaffold1_1 的模板序列字符数 6916724 - 模板序列行数1
=6916723
即p3_in.pl 把与misa文件ID序列名称对应的序列 全部作为模板连。就相当于p3in文件里的内容很大一部分为 基因组tsh_wgs.fna的碱基序列,单单scaffold1就出现了 3793次也即scaffold1查到到的ssr位点数。
6.3 PRIMER_PRODUCT_SIZE_RANGE=100-280
不知道这个是不是需要提前设置,还是在primer3相关的配置文件里配置这个内容
6.4 TARGET=2055,18
p3_in.pl 文件里默认把起始ssr的起始位点减去3 得到目标序列开始位点2055,ssr长度加 6 得到目标序列长度。
- 从基因组中提取只包含重复序列的一小段序列信息的bed文件。
7.1 cat tsh_wgs.fna.misa | awk -F "\t" 'NR >1{print$1"\t"$6-150"\t"$7+150}'> tsh_wgs_ssr.bed
在ssr序列起始终止位点向两端延伸150bp。
利用bed文件提取需要给基因组文件构建fai索引
7.2 samtools fadix tsh_wgs.fna
7.3 bedtools getfasta -fi tsh_wgs.fna -bed tsh_wgs_ssr.bed -fo tsh_wgs_ssr.fa
遇见问题Error: malformed BED entry at line 6643\. Start Coordinate detected that is < 0\. Exiting.
7.4
sed -n 6643p tsh_wgs_ssr.bed
scaffold1000 -30 281 #在scaffold 的较开始端有ssr位点,起始位置-150 变成了负数。#问题好多呀
7.5
sed -n 6644p tsh_wgs.fna.misa
scaffold1000 1 p2 (AG)6 12 120 131 #misa文件中ssr位点的起始位点为12
就是需要将bed文件中$2起始位置小于0的全部改为1,或在生成bed文件过程中将小于150的起始位置n-(n-1).
7.6 查到bed文件小于 起始位点数值小于0的行数
7.6.1 cat tsh_wgs_ssr.bed | awk ' $2 < 0{print $0}' |wc -l
78 #有78行起始位点小于0
7.6.2 cat tsh_wgs_ssr.bed | awk ' $2 < 0{print $1}' | wc -l
78 #只打印scaffold 这一列也是有78 行
7.6.3 cat tsh_wgs_ssr.bed | awk ' $2 < 0{print $1}' | sort | uniq -c |wc -l
77 #排除重复 只打印出哪些scaffold ,起始位点小于0,有77 条 扫scaffold,说明一条scaffold 上有2个位点的起始位置小于0
7.6.4 cat tsh_wgs_ssr.bed | awk ' $2 < 0{print $1}' | sort | uniq -c | awk '$1 >1{print$0}'
2 scaffold637 #scaffold637 里有2个ssr位点起始位置小于0
7.7 怎么样把小于0数值的改为1
cat tsh_wgs_ssr.bed |sed 's/-.*\t/1/g' | awk '$2 < 0 {print $0}' | wc -l
0 #匹配-负号内容,并改为1,没有查找到小于0的行,是因为把-.*\t替换为了1,即有的行只有2列
cat tsh_wgs_ssr.bed | sed 's/-.*\t/1\t/' | awk '$2<=1 {print($0)}' | wc -l
80 # 即起始位置小于等于1的行数有80行
grep -w -n "1" tsh_wgs_ssr.bed
23072:scaffold136 1 315
46722:scaffold2011 1 330 #起始位点等于1的有两行,说明可以通过匹配- 来更改小于0的起始位点。
grep -w -n "1" tsh_wgs_ssr.bed | wc -l
80 #更改后的等于1的有80 行
终于解决了,学了下正则表达式
- 借助bedtools提取bed文件的相应序列
8.1提取序列
bedtools getfasta -fi tsh_wgs.fna -bed tsh_wgs_ssr.bed - fo tsh_wgs_ssr.fa
Feature (scaffold1081:1618-1998) beyond the length of scaffold1081 size (1879 bp). Skipping.
Feature (scaffold1142:1594-1907) beyond the length of scaffold1142 size (1758 bp). Skipping.
.....................
#有的ssr位点在scaffold尾端,终止位置加150之后超过scaffold长度
8.2 处理终止位置大于序列总长度的信息行
复制 Feature 信息,将75行报错的信息内容更改,终止位置序列总bp数
$ cat *.bed > all.bed
#将 2个bed文件合并为一个bed文件。即原先的bed文件和修改后的bed文件。
8.3 提取all.bed文件的序列
$ bedtools getfasta -fi tsh_wgs.fna -bed all.bed -fo tsh_wgs_ssr.fa
$ grep -c ">" tsh_wgs_ssr.fa
127130
$ wc -l tsh_wgs.fna.misa
127131 tsh_wgs.fna.misa
tsh_wgs_ssr.fa 序列条数与 sh_wgs.fna.misa 的查找到的ID总数(减去表头)一样
8.4 使用misa鉴定位点
$ perl /home/Pomgroup/gdp/app/misa/misa.pl tsh_wgs_ssr.fa
16 Total number of sequences examined: 127130
17 Total size of examined sequences (bp): 40938166
18 Total number of identified SSRs: 155330
19 Number of SSR containing sequences: 127123
20 Number of sequences containing more than 1 SSR: 23222
21 Number of SSRs present in compound formation: 19525
奇怪,对提取的序列重新搜索SSRs与原基因组文件搜索到的SSRs数量不同
比基因组SSRs多出 155330 - 146445 = 8885个位点,
$ wc -L tsh_wgs.fna.misa
8383774 tsh_wgs.fna.misa 原misa文件中最长行,即复合卫星中间的插入大于100bp
8.4.1 看一下什么原因
@1
$ sed -n '58p' tsh_wgs.fna.misa
scaffold1 57 c (A)10gaaaaaaa(AAG)5gg(GA)10 55 99231 99285
#这个复合位点计算了插入片段片段的长度。
而以tsh_wgs.fna.misa第3388行为例
(base) [Pomgroup@localhost tshscaffold]$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$4}'|awk -F "[()]" '{print NF}'
5
(base) [Pomgroup@localhost tshscaffold]$ sed -n '3388p' tsh_wgs.fna.misa |awk -F "[()]" '{print NF}'
5
(base) [Pomgroup@localhost tshscaffold]$ sed -n '3388p' tsh_wgs.fna.misa |awk -F "[()]" '{print $1}'
scaffold1 3387 c
(base) [Pomgroup@localhost tshscaffold]$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$4}'|awk -F "[()]" '{print $1}'
$1字段开头有分隔符,所以可以直接以()为分隔符,查看某一列,
$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$4}'|awk -F "[()]" '{print$3}' |wc -m
617363 #可见插入片段的长度有617363 -1($2是7,为重复的次数)=617362 bp #有没有行尾符?
$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$5}'
34 #显示卫星的长度只有34bp,通过起始位点与终止位点的计算。
6299381 - 6299348 +1 =34
而不计算插入序列重复序列的总长度为2(CT)*7 + 1(T)*21=35 , 可能是把插入系列作为一个碱基?
$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$4}' |grep -o "[0-9]"
#查看重复序列的重复次数 数字
7
2
1
(CT)7....(T)21
$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$6}'
6299348 #起始位点
$ sed -n '3388p' tsh_wgs.fna.misa |awk '{print$7}'
6299381 #终止位点
misa的坐标系统为“1-start, fully-closed” system,即从重复开始的那个碱基计算位点位置到终止碱基位点。
@2 查看3388行所在scaffold1的bp数
grep -n "scaffold" tsh_wgs.fna |less -SN
1 1:>scaffold1
2 86462:>scaffold10
即碱基位于2-86461行
$ sed -n '2,86461p' tsh_wgs.fna |wc -l
86460
$ sed -n '2,86461p' tsh_wgs.fna |wc -m
7003183
7003183-86460=6919723 scaffold1长度
??没有找到插入序列这么长的原因。
- 处理misa文件
$ perl p3_in.pl tsh_wgs_ssr.fa.misa
135805 records created.
$ wc -l tsh_wgs_ssr.fa.misa
135806 tsh_wgs_ssr.fa.misa
#所有行都准确识别
10 .使用primer3 批量设计引物
前段时间学的已经忘了差不多啦,反正是有各种参数,
感谢作者对各种参数的解释
primer3 的使用参数
/home/Pomgroup/gdp/app/primer/primer3-2.4.0/src/primer3_core [--format_output] [--default_version=1|--default_version=2] [--io_version=4] [--p3_settings_file=<file_path>] [--echo_settings_file] [--strict_tags] [--output=<file_path>] [--error=<file_path>] [input_file]
$ /home/Pomgroup/gdp/app/primer/primer3-2.4.0/src/primer3_core --p3_settings_file=/home/Pomgroup/gdp/app/primer/primer3-2.4.0/settings_files/p3_th_settings.txt --output= tsh_ssr.p3out --error= tsh.error tsh_wgs_ssr.fa.p3in
#程序没有反应,可能是因为--output=后于输出文件之间有空格。更改后运行正常
10.1 p3_out.pl 处理设计好的引物
perl tsh_ssr.p3out tsh_wgs_ssr.fa.misa
Semicolon seems to be missing at tsh_ssr.p3out line 1.
Can't modify constant item in scalar assignment at tsh_ssr.p3out line 1, near "s caffold1:"
syntax error at tsh_ssr.p3out line 1, near "scaffold1:"
Identifier too long at tsh_ssr.p3out line 2.
#perl后没有 p3_out.pl
$ perl p3_out.pl tsh_ssr.p3out tsh_wgs_ssr.fa.misa
Primer modelling was successful for sequences.
Primer modelling failed for 135805 sequences.
# 产生 的tsh_ssr.result文件为空
用师兄的python程序可以将文件整理到一个结果
$ wc -l tsh_ssr.result
97156 tsh_ssr.result #相当于结果的总共的位点数有97155个
而
$ grep -c "SEQUENCE_ID" tsh_ssr.p3out
135805 #文件有135805条模板序列
同时
$ wc -l tsh_wgs_ssr.fa.misa
135806 tsh_wgs_ssr.fa.misa #除去表头misa文件也有135805行
135805 - 97155 = 35650 个位点没有合并到结果文件
可见misa在scaffold1011查到到一个位点,而result文件里没有出现scaffold1011里的位点,查看p3out文件里,scaffold1011引物设计情况
scaffold1011查到到的唯一没有设计出引物
SEQUENCE_ID=scaffold1011:1563-1892_1
SEQUENCE_TEMPLATE=AATAATACTGAAAGCTTGTGGTTAGATATGCTGATCACTTGTCCAAATCTAGAGGACTACtaattgatattatattaCCGAGCTCGTGGTTAGATGGATACATGGACAATTTTATGTAGAAGACTACTAATTGATTTATTACtgatgtgatatatatatatatatatatatatatatatctattacaCCAAGCGCGATGTACTCAAATCTTATGATATTTTACCCGAAAGGATTCAAGGGAGAAAACCCAGATCTATCTGTCTACATATGTTCGATATTACAGATACTCTTTTTTGATGAATGTGTATTATTGATTCTCATTATCATTA
PRIMER_PRODUCT_SIZE_RANGE=100-280
TARGET=147,36
PRIMER_MAX_END_STABILITY=250
PRIMER_LEFT_EXPLAIN=considered 1380, GC content failed 998, low tm 310, high tm 23, high hairpin stability 13, ok 36
PRIMER_RIGHT_EXPLAIN=considered 1380, GC content failed 1080, low tm 243, high tm 15, ok 42
PRIMER_PAIR_EXPLAIN=considered 1512, unacceptable product size 1512, ok 0
PRIMER_LEFT_NUM_RETURNED=0
PRIMER_RIGHT_NUM_RETURNED=0
PRIMER_INTERNAL_NUM_RETURNED=0
PRIMER_PAIR_NUM_RETURNED=0
=
查看scaffold110 在misa文件与result文件前后变化
#将 tsh_ssr.p3out 文件中 SEQUENCE_ID= 中有scaffold110的内容提取的文件
$ sed -n '/scaffold110:/p' tsh_ssr.p3out |awk -F "=" '{print$2}' > p3out_scaffold110_ssr.txt
#将result文件中 第一列ID为 scaffold110 提取出来
$ sed -n '/scaffold110:/p' tsh_ssr.result |awk -F "\t" '{print$1}' > result_scaffold110_ssr.txt
#查看两个文件中重复与没有重复的内容,统计到的个数为1的就是id就是在最后整理到的时候没有整理到最后的result文件里id,可是是这些id没有设计出引物。
$ cat result_scaffold110_ssr.txt p3out_scaffold110_ssr.txt |sort|uniq -c
#将两个文件没有重复的内容提取到 p3out_result_miss.txt
$ cat result_scaffold110_ssr.txt p3out_scaffold110_ssr.txt |sort|uniq -c | awk '$1~/1/{print$2}' > p3out_result_miss.txt
#统计下miss掉的id个数
wc -l p3out_result_miss.txt
71 p3out_result_miss.txt #等于上边scaffold110 在两个文件中的差值即316-245=71
#查看这71个id因为什么原因没有被统计到结果文件中,根据p3out_result_miss.txt里的id名称在p3out文件里提取对应的所有内容。#不会只能现学了,用python?#(大写的搁置)
统计result文件内 各scaffold的包括位点数
$cat tsh_ssr.result | awk -F ":" 'NR>1{print $1}' | sort | uniq -c > result.statistics
以tsh_wgs_ssr.fa.misa 文件中scaffold1000为例只检测到3个位点?
$ cat tsh_wgs_ssr.fa.misa | awk -F ":" '$1~"scaffold1000" {print$0}'
scaffold1000:1-281 1 p2 (AG)6 12 119 130
scaffold1000:131-480 1 c (CTC)8cttcctccatt(TTC)5 50 150 199
scaffold1000:1414-1733 1 p2 (TC)10 20 150 169
$ grep -c "scaffold1000:" tsh_wgs_ssr.fa #序列文件中scaffold1000也是只有3个位点
$ grep -c "scaffold1000" tsh_wgs.fna.misa
3 #也是3个位点,是因为scaffold1000序列短,一开始以为scaffold都很长呢啊啊啊啊,
10 使用e-PCR检测引物的有效性
$ nohup /home/Pomgroup/gdp/app/epcr/e-PCR-2.3.12/e-PCR /home/Pomgroup/gdp/ssr/tsh/tshscaffold/ssrfa/ssrepcr/tsh_ssr_result_epcr.txt D=100-500 /home/Pomgroup/gdp/ssr/tsh/tshscaffold/tsh_wgs.fna N=2 G=2 T=3 > tsh_ssr_result_afepcr.txt &
结果文件内容为
scaffold1 scaffold1:5934613-5934928_1 + 5934661 5934761 101/100-500 0 0
scaffold1 scaffold1:2910453-2910765_1 + 2910543 2910724 182/100-500 0 0
scaffold1 scaffold1:879458-879772_1 + 879469 879675 207/100-500 0 0
#加号一列有个别减号,还不清楚+ - 的含义,可能是代表是否扩增成功。
e-PCR的参数也不懂,先学一下e-PCR 的说明书。
10.1 e-PCR的参数
M=## Margin (default 50)
W=## Wordsize (default 7)
N=## Max mismatches allowed (default 0) 允许的最大错配,(是碱基错
配么不清楚)
G=## Max indels allowed (default 0) 允许的indel最大错配
F=## Use ## discontinuos words
O=## Set output file to ##
T=## Set output format (1..4) 输出的文件格式
D=##-## Set default sts size #sts长度
P=+- Postprocess hits on/off
V=+- Verbose on/Off
A=a|f Use presize alignmens (only if gaps>0), slow
a - Allways or f - as Fallback
X=+- Use 5'-end lowercase masking of primers (default -)
U=+- Uppercase all primers (default -)
-mid Same as T=2
compat-options 可以出现在命令行的任何地方
#EPCR版本
$ /home/Pomgroup/gdp/app/epcr/e-PCR-2.3.12/e-PCR -V
e-PCR cmdline tool version 2.3.12
即epcr将result文件中的left 引物与基因组比对,将result文件中的right引物的反向互补序列与基因组比对,如果完全比对则表明可以扩增。
对提取后的
11 筛查唯一结果
X 先挑选设计出的类型为 重复类型P3 的引物,序列较长的位点设计引物,看设计出来引物是否可以扩增出条带。
查看重复序列为p3的位点
$ cat tsh_ssr.result |awk 'BEGIN{FS=OFS="\t"}$3~/p3/{print$0}' | awk 'BEGIN{FS=OFS="\t"}$5>15{print$0}' |less -SN
挑选不同scaffold上的位点总共20对,
将序列全部改为大写重新设计引物,一是epcr对引物的大小写有要求,二是第一次的引物设计参数设置不合理吧,需要重新设计引物
采用默认参数设置
$ nohup /home/Pomgroup/gdp/app/primer/primer3-2.4.0/src/primer3_core --default_version=2 --output=tsh_wgs_ssr_upper.fa.p3out --error= tsh_wgs_ssr_upper.fa.error tsh_wgs_ssr_upper.fa.p3in && echo "well done" || echo "failure" &
挑选p3类型引物20对选择进行预实验
$ cat tsh_wgs_ssr_upper.result |awk 'BEGIN{FS=OFS="\t"}$3~/p3/{print$0}'|awk 'BEGIN{FS=OFS="\t"}$5>15{print$0}' | > tsh_wgs_ssr_upper.p3great15.result
bedtools getfasta -fi tns.fna -bed tns_misa_ssr_big0.bed -fo tns_misa_ssr.fna
删除起始位置小于0的 bed文件行
关于引物设计
上下为反向互补序列,上边为5-3方向,下边为3-5方向
中间标记颜色的为SSR重复序列,引物设计是根据位点两端序列进行设计,并不是从位点前的第一个碱基往左数多少个碱基作为引物,因为引物设计有要求,所有引物最后一个碱基往前移动,并调整引物长度,计算后判断哪段序列作为引物合适。所有ssr引物扩增的产物长度是引物长度,+ssr序列长度,+引物3端最后一个碱基据ssr序列5端的长度
misa生成的后缀为misa的文件时每个位点的信息,statistics文件是一个总概括文件,statistics文件里的统计的1-6重复型ssr的数量是把复合微卫星中的单个微卫星计算进去,加入在misa后缀文件里有p1(单核苷酸重复)类型经手动查找计数共10个,在statistics文件里p1(单核苷酸重复)类型有15个,说明所有复合微卫星内含有5个单独的单核苷酸重复ssr。