植物长链非编码RNA（lncRNA）鉴定实例（拟南芥数据）

SEQLIBS=(SRR8428909 SRR8428908 SRR8428906 SRR8428905 )

for seqlib in ${SEQLIBS[@]}; do
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR842/00${seqlib:9:10}/${seqlib}/${seqlib}_1.fastq.gz
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR842/00${seqlib:9:10}/${seqlib}/${seqlib}_2.fastq.gz
done

bash download_raw_data.sh

bgzip -d SRR8428909_1.fastq.gz
bgzip -d SRR8428909_2.fastq.gz

mv SRR8428909_1.fastq wt_Rep1_R1.fastq
mv SRR8428909_2.fastq wt_Rep1_R2.fastq

bgzip -d SRR8428908_1.fastq.gz
bgzip -d SRR8428908_2.fastq.gz

mv SRR8428908_1.fastq wt_Rep2_R1.fastq
mv SRR8428908_2.fastq wt_Rep2_R2.fastq

bgzip -d SRR8428906_1.fastq.gz
bgzip -d SRR8428906_2.fastq.gz

mv SRR8428906_1.fastq EE_Rep1_R1.fastq
mv SRR8428906_2.fastq EE_Rep1_R2.fastq

bgzip -d SRR8428905_1.fastq.gz
bgzip -d SRR8428905_2.fastq.gz

mv SRR8428905_1.fastq EE_Rep2_R1.fastq
mv SRR8428905_2.fastq EE_Rep2_R2.fastq

SEQLIBS=(EE_Rep1 EE_Rep2  wt_Rep1 wt_Rep2 )

for seqlib in ${SEQLIBS[@]}; do
    fastp -i ${seqlib}_R1.fastq -I ${seqlib}_R2.fastq -o ${seqlib}_clean_R1.fastq -O ${seqlib}_clean_R2.fastq
done

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
bgzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
mv  Arabidopsis_thaliana.TAIR10.dna.toplevel.fa At.fa
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.40.gff3.gz
bgzip -d Arabidopsis_thaliana.TAIR10.40.gff3.gz 
mv Arabidopsis_thaliana.TAIR10.40.gff3 At.gff3

mkdir reference
mv At* reference
cd reference
bowtie2-build At.fa At

cd ../
tophat2 -p 8 -I 5000 -G reference/At.gff3 -o wt1_thout reference/At ../wt_Rep1_clean_R1.fastq ../wt_Rep1_clean_R2.fastq
tophat2 -p 12 -I 5000 -G reference/At.gff3 -o wt2_thout reference/At ../wt_Rep2_clean_R1.fastq ../wt_Rep2_clean_R2.fastq

tophat2 -p 8 -I 5000 -G reference/At.gff3 -o EE1_thout reference/At ../EE_Rep1_clean_R1.fastq ../EE_Rep1_clean_R2.fastq
tophat2 -p 8 -I 5000 -G reference/At.gff3 -o EE2_thout reference/At ../EE_Rep2_clean_R1.fastq ../EE_Rep2_clean_R2.fastq

cufflinks -p 4 -g reference/At.gff3 -I 5000 -o wt1_clout wt1_thout/accepted_hits.bam
cufflinks -p 4 -g reference/At.gff3 -I 5000 -o wt2_clout wt2_thout/accepted_hits.bam
cufflinks -p 4 -g reference/At.gff3 -I 5000 -o EE1_clout EE1_thout/accepted_hits.bam
cufflinks -p 4 -g reference/At.gff3 -I 5000 -o EE2_clout EE2_thout/accepted_hits.bam

find . -name transcripts.gtf > assemblies.txt
cuffmerge -p 4 -g reference/At.gff3 -s reference/At.fa assemblies.txt

cuffdiff -o diff_out -b reference/At.fa -p 8 -T -L EE,wt -u merged_asm/merged.gtf EE1_thout/accepted_hits.bam,EE2_thout/accepted_hits.bam wt1_thout/accepted_hits.bam,wt2_thout/accepted_hits.bam

cat merged_asm/merged.gtf | grep 'class_code "[uiox]"' > selected.gtf
gffread -w selected.fa -g reference/At.fa selected.gtf
cat selected.fa | grep ">" | wc -l

import sys
from Bio import SeqIO

inputfasta = sys.argv[1]
min_len = int(sys.argv[2])
outputfasta = sys.argv[3]

i = 0
j = 0

fw = open(outputfasta,'w')

for rec in SeqIO.parse(inputfasta,'fasta'):
	i += 1
	if len(rec.seq) > min_len:
		j += 1
		fw.write(">%s\n%s\n"%(rec.id,str(rec.seq)))
		
print("The total number of sequence is: ",i)
print("Number of sequence retained is: ",j)

python .\select_seq_according_to_seq_length.py .\selected.fa 200 output.fasta

cat result_cpc2.txt | grep 'noncoding' | cut -f 1 > non-coding-transcript.txt
wc -l non-coding-transcript.txt