降解组测序原理： AGO蛋白会在miRNA与mRNA互补区域的第10位碱基（从miRNA5’端开始数）处切割靶基因mRNA序列，靶基因序列分为5’片段和3’片段，其中3’片段由于含有游离的5’磷酸基团，会被RNA酶连接上5’Adaptor，5’片段和未被剪切的mRNA由于5端含有帽子结构无法被加上接头而不会进行测序，5’Adaptor序列中含有一种内切酶MmeI的识别位点，这个酶切割的位点是在识别结合位点的20-30bp后，切割合成的双链cDNA后加上3’Adaptor进行测序；降解组序列长度是50 nt左右，故测序得到的序列上有一段 3’接头序列。
联川新的建库流程好像不需要酶切（http://www.lc-bio.com/Techservices/show-133.html）
用降解组数据预测phasiRNAs靶基因需要准备以下3个文件：
1、去除接头的降解组fasta序列文件（ID简短，不要有空格，文件名也不要有空格）；
2、小RNA query文件（ID简短，不要有空格）；
3、转录本序列作为靶基因库；
这3个文件的获取方式：
1、NCBI下载降解组数据文件；
将SRA格式的文件转换成fastq格式，以小RNA去接头的方法去除降解组数据的接头；降解组数据不需要去冗余；
2、根据课题需求获取含有miRNA靶位点的PHAS locus序列（靶位点在头或者尾），miRNA trigger PHAS locus产生phsiRNAs，其切割位点一般是靶位点的第10位碱基（miRNA的5’端开始计数）；写脚本取出以21nt相位切割的方式产生的siRNAs，并筛选保留在sRNA数据中有表达丰度的siRNAs;
#脚本一：获取phasiRNAs序列的脚本 get_phasiRNAs_seq.py
#!/usr/bin/env python
# -*- coding:utf-8 -*-
seq = input('输入PHAS序列（靶位点所在序列）:')
f = open('phasiRNA.seq','w')
complement = {'C': 'G', 'G': 'C', 'T': 'A', 'A': 'T'}
com_seq = ''
for i in list(seq):
com_seq += complement[i]
rev_com_seq = com_seq[::-1]
rc_seq = rev_com_seq
print(rc_seq,)
start1 = 12
start2 = 11
START1 = -15
START2 = -14
for i in range(20):
end1 = start1 + 21
phasiRNA1 = seq[start1:end1]
start1 = end1
print(phasiRNA1,i,'10th',file = f)
end2 = start2 + 21
phasiRNA2 = seq[start2:end2]
start2 = end2
print(phasiRNA2,i,'11th',file = f)
END1 = START1 - 21
siRNA1 = rc_seq[END1:START1]
START1 = END1
print(siRNA1,'-%d'%(i),'10th',file = f)
END2 = START2 - 21
siRNA2 = rc_seq[END2:START2]
START2 = END2
print(siRNA2,'-%d'%(i),'11th',file = f)
f.close()
python get_phasiRNAs_seq.py
#输入靶位点靶向的链序列
#脚本二：筛选在sRNA数据中有表达的phasiRNAs并以fasta格式输出 extract_abund.py
#!/usr/bin/env python
# -*- coding:utf-8 -*-
import sys
File1 = sys.argv[1]
File2 = sys.argv[2]
collapsed_fasta = open(File1,'r')
siRNA_seq = open(File2,'r')
prefix1 = sys.argv[1].split("_")[0]
prefix2 = sys.argv[2].split(".")[0]
output_file = open ('%s_%s.seq'%(prefix1,prefix2), 'w')
Dict = dict()
i = 0
for line1 in collapsed_fasta:
line1 = line1.rstrip()
if not line1: break
if '>' in line1:
abund = line1.split('-')[1]
else:
Dict[line1] = abund
for line2 in siRNA_seq:
if not len(line2):continue
line2 = line2.rstrip()
list2 = line2.split()
for key, value in Dict.items():
if key == list2[0]:
i += 1
print(">%d-%s-%s-%s"%(i,value,list2[1],list2[2]),file = output_file)
print(key,file = output_file)
collapsed_fasta.close()
siRNA_seq.close()
output_file.close()
python extract_abund.py SRR4010495_mc.fa PtncRNAa.seq
3、Phytozome下载物种的转录组数据
用TBtools简化转录组序列ID-Fasta ID Simplifier
以上3个文件准备完成后，用CleaveLand进行靶基因预测
CleaveLand4.pl -e SRR4010497.fastq.trimmed -u PtncRNAa_phasiRNAs.fasta -n Ptrichocarpa_444_v3.1.cds_primaryTranscriptOnly.fa.sim -o plot/ > PtncRNAa_SRR4010497_PARE.txt
#产生的文件如下：
#SRR4010497.fastq.trimmed_dd.txt
#PtncRNAa_phasiRNAs.fasta_GSTAr.txt
#PtncRNAa_SRR4010497_PARE.txt
#Plot（含有预测T-plot图的pdf文件）
其中PtncRNAa_SRR4010497_PARE.txt文件含有一个siRNA对应的一个靶基因的详细信息，较为重要的信息有Alignment score和Category；
Alignment score计算原则： mismatch或bulge罚分为1， G-U摇摆罚分0.5，在小RNA 5‘端开始数第2-13碱基中所有罚分加倍；AS越小越可信。
Category的分类：
Category 4: 只有1个read的位置
Category 3: >1 read,但是低于或等于覆盖在转录本上的reads的平均深度
Category 2: >1 read,高于覆盖在转录本上的reads的平均深度，但不是最大的
Cateogry 1: >1 read,等于覆盖在转录本上的reads的最大深度，且最大值的位置不止1个
Cateogry 0: >1 read, 等于覆盖在转录本上的reads的最大深度，且最大值的位置只有1个
Category越小越可信。
以下为一个较为可信的靶位点：siRNA 14-43靶向切割转录本Potri.005G227300.1的第105位碱基，红点为切割位点。
用IGV查看T-plot中的每个切割位点是否有降解组reads
STAR --runThreadN 30 --runMode alignReads --genomeDir ./genomeIndex --readFilesIn SRR4010497.fastq.trimmed --outFileNamePrefix SRR4010497_ --alignIntronMax 20000 --alignMatesGapMax 25000 --outFilterMultimapNmax 50 --outSAMtype BAM SortedByCoordinate
samtools index SRR4010497_Aligned.sortedByCoord.out.bam
IGV查看切割产生reads的位置与T-plot的位置信息相对应，可以作为验证
提取该转录本的序列用Swissprot进行注释
CleaveLand使用注意事项：
https://github.com/MikeAxtell/CleaveLand4
1、转录组序列的ID行不要太长，不要有空格，转录组序列文件的文件名也不要有空格；
e.g. ">AT1G12345" is good, ">AT1G12345 | this is my favorite gene | it is awesome" is not.
2、降解组序列文件要去接头，去接头方式与小RNA数据一样，且降解组数据不需要去冗余，但听说联川现在的测序数据不用去接头了？
3、小RNA序列的ID行同样需要简短不含空格，DNA序列或者RNA序列都可以；
e.g. ">ath-miR169a" is good, ">ath-miR169a MIMAT0000200 Arabidopsis thaliana miR169a" is not.
4、如果在miRBase中下载了一个物种中同源miRNA基因的完全一致的成熟序列，先去除冗余；
5、不要用全基因组数据，用CDS转录本数据（考虑到假阳性及内存问题）；
6、Query和转录本序列都不要含有简并碱基，含有简并碱基的转录本会被忽略；
7、Query的长度为15-26 nt；
8、CleaveLand4只会寻找比对小RNA第10位碱基切割的证据；目前还没有直接的证据能证明AGO蛋白能切割小RNA除第10位之外的位置。
